{ "cells": [ { "cell_type": "markdown", "id": "4ee595c2", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Monte Carlo and Root Finding" ] }, { "cell_type": "markdown", "id": "d0028b39", "metadata": { "slideshow": { "slide_type": "-" }, "tags": [ "remove-cell" ] }, "source": [ "**CS1302 Introduction to Computer Programming**\n", "___" ] }, { "cell_type": "markdown", "id": "711e203e", "metadata": { "slideshow": { "slide_type": "subslide" }, "tags": [] }, "source": [ "## The Monty-Hall Game" ] }, { "cell_type": "code", "execution_count": 1, "id": "979807da", "metadata": { "slideshow": { "slide_type": "-" }, "tags": [ "hide-input" ] }, "outputs": [ { "data": { "text/html": [ "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%%html\n", "" ] }, { "cell_type": "markdown", "id": "eb2223c1", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "**Is it better to change the initial pick? What is the chance of winning if we change?**" ] }, { "cell_type": "markdown", "id": "7221934c", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "*Hint:* There are two doors to choose from, and only one of the doors has treasure behind." ] }, { "cell_type": "markdown", "id": "d2da94c3", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "Let's use the following program to play the game a couple of times." ] }, { "cell_type": "code", "execution_count": 2, "id": "35dca6ed", "metadata": { "code_folding": [], "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "import random\n", "\n", "\n", "def play_monty_hall(num_of_doors=3):\n", " \"\"\"Main function to run the Monty Hall game.\"\"\"\n", " doors = {str(i) for i in range(num_of_doors)}\n", " door_with_treasure = random.sample(doors, 1)[0]\n", "\n", " # Input initial pick of the door.\n", " while True:\n", " initial_pick = input(f'Pick a door from {\", \".join(sorted(doors))}: ')\n", " if initial_pick in doors:\n", " break\n", "\n", " # Open all but one other door. Opened door must have nothing.\n", " doors_to_open = doors - {initial_pick, door_with_treasure}\n", " other_door = (\n", " door_with_treasure\n", " if initial_pick != door_with_treasure\n", " else doors_to_open.pop()\n", " )\n", " print(\"Door(s) with nothing behind:\", *doors_to_open)\n", "\n", " # Allow player to change the initial pick the other (unopened) door.\n", " change_pick = (\n", " input(f\"Would you like to change your choice to {other_door}? [y/N] \").lower()\n", " == \"y\"\n", " )\n", "\n", " # Check and record winning.\n", " if not change_pick:\n", " mh_stats[\"# no change\"] += 1\n", " if door_with_treasure == initial_pick:\n", " mh_stats[\"# win without changing\"] += 1\n", " return print(\"You won!\")\n", " else:\n", " mh_stats[\"# change\"] += 1\n", " if door_with_treasure == other_door:\n", " mh_stats[\"# win by changing\"] += 1\n", " return print(\"You won!\")\n", " print(f\"You lost. The door with treasure is {door_with_treasure}.\")\n", "\n", "\n", "mh_stats = dict.fromkeys(\n", " (\"# win by changing\", \"# win without changing\", \"# change\", \"# no change\"), 0\n", ")\n", "\n", "\n", "def monty_hall_statistics():\n", " \"\"\"Print the statistics of the Monty Hall games.\"\"\"\n", " print(\"-\" * 30, \"Statistics\", \"-\" * 30)\n", " if mh_stats[\"# change\"]:\n", " print(\n", " f\"% win by changing: \\\n", " {mh_stats['# win by changing'] / mh_stats['# change']:.0%}\"\n", " )\n", " if mh_stats[\"# no change\"]:\n", " print(\n", " f\"% win without changing: \\\n", " {mh_stats['# win without changing']/mh_stats['# no change']:.0%}\"\n", " )" ] }, { "cell_type": "code", "execution_count": 3, "id": "90457207", "metadata": { "slideshow": { "slide_type": "subslide" }, "tags": [ "remove-output" ] }, "outputs": [ { "ename": "StdinNotImplementedError", "evalue": "raw_input was called, but this frontend does not support input requests.", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mStdinNotImplementedError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mplay_monty_hall\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2\u001b[0m \u001b[0mmonty_hall_statistics\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m\u001b[0m in \u001b[0;36mplay_monty_hall\u001b[0;34m(num_of_doors)\u001b[0m\n\u001b[1;32m 9\u001b[0m \u001b[0;31m# Input initial pick of the door.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 10\u001b[0m \u001b[0;32mwhile\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 11\u001b[0;31m \u001b[0minitial_pick\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0minput\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34mf'Pick a door from {\", \".join(sorted(doors))}: '\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 12\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0minitial_pick\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mdoors\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 13\u001b[0m \u001b[0;32mbreak\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m~/my-conda-envs/jb/lib/python3.8/site-packages/ipykernel/kernelbase.py\u001b[0m in \u001b[0;36mraw_input\u001b[0;34m(self, prompt)\u001b[0m\n\u001b[1;32m 843\u001b[0m \"\"\"\n\u001b[1;32m 844\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_allow_stdin\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 845\u001b[0;31m raise StdinNotImplementedError(\n\u001b[0m\u001b[1;32m 846\u001b[0m \u001b[0;34m\"raw_input was called, but this frontend does not support input requests.\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 847\u001b[0m )\n", "\u001b[0;31mStdinNotImplementedError\u001b[0m: raw_input was called, but this frontend does not support input requests." ] } ], "source": [ "play_monty_hall()\n", "monty_hall_statistics()" ] }, { "cell_type": "markdown", "id": "03b7905d", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "You may also [play the game online](https://math.ucsd.edu/~crypto/Monty/monty.html)." ] }, { "cell_type": "markdown", "id": "8a13e7c7", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "To get a good estimate of the chance of winning, we need to play the game many times. \n", "We can write a Monty-Carlo simulation instead." ] }, { "cell_type": "code", "execution_count": 4, "id": "004ae038", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "code", "checksum": "e8cf784178caa46c7fbd46c61ba83e54", "grade": false, "grade_id": "monty-hall", "locked": true, "schema_version": 3, "solution": false, "task": false }, "slideshow": { "slide_type": "-" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Door with treasure: 1 3 1 3 1 2 3 3 1 1 ...\n", " Initial pick: 1 2 1 2 1 1 1 1 2 2 ...\n" ] } ], "source": [ "# Do not change any variables defined here, or some of the tests may fail.\n", "import numpy as np\n", "\n", "np.random.randint?\n", "\n", "np.random.seed(0) # for reproducible result\n", "num_of_games = int(10e7)\n", "door_with_treasure = np.random.randint(1, 4, num_of_games, dtype=np.uint8)\n", "initial_pick = np.random.randint(1, 4, num_of_games, dtype=np.uint8)\n", "\n", "print(f\"{'Door with treasure:':>19}\", *door_with_treasure[:10], \"...\")\n", "print(f\"{'Initial pick:':>19}\", *initial_pick[:10], \"...\")" ] }, { "cell_type": "markdown", "id": "a940b7ae", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- `door_with_treasure` stores as 8-bit unsigned integers `uint8` the door numbers randomly chosen from $\\{1, 2, 3\\}$ as the doors with treasure behind for a number `num_of_games` of Monty-Hall games.\n", "- `initial_pick` stores the initial choices for the different games." ] }, { "cell_type": "markdown", "id": "8f424bc3", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "If players do not change their initial pick, the chance of winning can be estimated as follows:" ] }, { "cell_type": "code", "execution_count": 5, "id": "3ae296cd", "metadata": { "code_folding": [] }, "outputs": [ { "data": { "text/plain": [ "0.333291" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def estimate_chance_of_winning_without_change(door_with_treasure, initial_pick):\n", " \"\"\"Estimate the chance of winning the Monty Hall game without changing\n", " the initial pick using the Monte Carlo simulation of door_with_treasure\n", " and initial_pick.\"\"\"\n", " count_of_win = 0\n", " for x, y in zip(door_with_treasure, initial_pick):\n", " if x == y:\n", " count_of_win += 1\n", " return count_of_win / n\n", "\n", "\n", "n = num_of_games // 100\n", "estimate_chance_of_winning_without_change(door_with_treasure[:n], initial_pick[:n])" ] }, { "cell_type": "markdown", "id": "af3077a4", "metadata": {}, "source": [ "However, the above code is inefficient and takes a long time to run. You may try running it on the entire sequences of `door_with_treasure` and `initial_pick` but **DO NOT** put the code in your notebook, as the server refuses to autograde notebooks that take too much time or memory to run." ] }, { "cell_type": "markdown", "id": "15f7d381", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "A simpler and also more efficient solution with well over 100 times speed up is as follows:" ] }, { "cell_type": "code", "execution_count": 6, "id": "0b22c164", "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/plain": [ "0.33332177" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def estimate_chance_of_winning_without_change(door_with_treasure, initial_pick):\n", " \"\"\"Estimate the chance of winning the Monty Hall game without changing\n", " the initial pick using the Monte Carlo simulation of door_with_treasure\n", " and initial_pick.\"\"\"\n", " return (door_with_treasure == initial_pick).mean()\n", "\n", "\n", "estimate_chance_of_winning_without_change(door_with_treasure, initial_pick)" ] }, { "cell_type": "markdown", "id": "68dcf80e", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "The code uses the method `mean` of `ndarray` that computes the mean of the `numpy` array. \n", "In computing the mean, `True` and `False` are regarded as `1` and `0` respectively, as illustrated below." ] }, { "cell_type": "code", "execution_count": 7, "id": "49398876", "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "True + True == 2\n", "True + False == 1\n", "False + True == 1\n", "False + False == 0\n" ] } ], "source": [ "for i in True, False:\n", " for j in True, False:\n", " print(f\"{i} + {j} == {i + j}\")" ] }, { "cell_type": "markdown", "id": "4e5d339e", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "**Exercise** Define the function `estimate_chance_of_winning_by_change` same as `estimate_chance_of_winning_without_change` above but returns the estimate of the chance of winning by changing the initial choice instead. Again, *implement efficiently or jupyterhub may refuse to autograde your entire notebook*.\n", "\n", "*Hint:* Since there are only two unopened doors at the end of each game, a player will win by changing the initial pick if the initially picked door is not the door with treasure behind." ] }, { "cell_type": "code", "execution_count": 8, "id": "6593e5ba", "metadata": { "deletable": false, "nbgrader": { "cell_type": "code", "checksum": "6b9be14256e67db2acf4f2ae073a57ae", "grade": false, "grade_id": "switch", "locked": false, "schema_version": 3, "solution": true, "task": false }, "slideshow": { "slide_type": "-" }, "tags": [ "remove-output" ] }, "outputs": [], "source": [ "def estimate_chance_of_winning_by_change(door_with_treasure, initial_pick):\n", " \"\"\"Estimate the chance of winning the Monty Hall game by changing\n", " the initial pick using the Monte Carlo simulation of door_with_treasure\n", " and initial_pick.\"\"\"\n", " # YOUR CODE HERE\n", " raise NotImplementedError()" ] }, { "cell_type": "code", "execution_count": 9, "id": "38462e40", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "code", "checksum": "29d0c3318c2c66dbf7c1792859a908ca", "grade": true, "grade_id": "test-switch", "locked": true, "points": 1, "schema_version": 3, "solution": false, "task": false }, "slideshow": { "slide_type": "-" }, "tags": [ "hide-input", "remove-output" ] }, "outputs": [ { "ename": "NotImplementedError", "evalue": "", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mNotImplementedError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0;31m# tests\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2\u001b[0m assert np.isclose(\n\u001b[0;32m----> 3\u001b[0;31m \u001b[0mestimate_chance_of_winning_by_change\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdoor_with_treasure\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;36m10\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0minitial_pick\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;36m10\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 4\u001b[0m \u001b[0;36m0.7\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 5\u001b[0m )\n", "\u001b[0;32m\u001b[0m in \u001b[0;36mestimate_chance_of_winning_by_change\u001b[0;34m(door_with_treasure, initial_pick)\u001b[0m\n\u001b[1;32m 4\u001b[0m and initial_pick.\"\"\"\n\u001b[1;32m 5\u001b[0m \u001b[0;31m# YOUR CODE HERE\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 6\u001b[0;31m \u001b[0;32mraise\u001b[0m \u001b[0mNotImplementedError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", "\u001b[0;31mNotImplementedError\u001b[0m: " ] } ], "source": [ "# tests\n", "assert np.isclose(\n", " estimate_chance_of_winning_by_change(door_with_treasure[:10], initial_pick[:10]),\n", " 0.7,\n", ")" ] }, { "cell_type": "code", "execution_count": 10, "id": "a1922f4e", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "code", "checksum": "4d5596f079ce618f7b208bfaa5df26e5", "grade": true, "grade_id": "h_test-switch", "locked": true, "points": 1, "schema_version": 3, "solution": false, "task": false }, "tags": [ "remove-cell" ] }, "outputs": [], "source": [ "# hidden tests" ] }, { "cell_type": "markdown", "id": "ad57bbb0", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "## Solving a 3-by-3 system of linear equations" ] }, { "cell_type": "markdown", "id": "887a7762", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "`numpy` has a module `linalg` for linear algebra, and the module provides a function called `solve` that can solve a system of linear equations. For the example in the lecture\n", "\n", "$$\n", "\\begin{aligned}\n", "2 x_0 + 2 x_1 &= 1\\\\\n", "2 x_1 &= 1,\n", "\\end{aligned}\n", "$$\n", "we can obtain the solution as follows:" ] }, { "cell_type": "code", "execution_count": 11, "id": "e06a58f6", "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "np.linalg.solve?\n", "A = np.array([[2.0, 2], [0, 2]])\n", "b = np.array([1.0, 1])\n", "x = np.linalg.solve(A, b)" ] }, { "cell_type": "markdown", "id": "9e4cf41e", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "As explained in the lecture, the arguments `A` and `b` are obtained from the matrix form of the system of linear equations:\n", "\n", "$$\n", "\\underbrace{\n", "\\begin{bmatrix}\n", "2 & 2\\\\\n", "0 & 2\n", "\\end{bmatrix}}_{\\mathbf{A}}\n", "\\underbrace{\n", "\\begin{bmatrix}\n", "x_0\\\\ x_1\n", "\\end{bmatrix}}_{\\mathbf{x}}\n", "= \n", "\\underbrace{\n", "\\begin{bmatrix}\n", "1 \\\\ 1\n", "\\end{bmatrix}\n", "}_{\\mathbf{b}}\n", "$$" ] }, { "cell_type": "markdown", "id": "ce5eca1f", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "However, the function returns an error when there is no unique solution." ] }, { "cell_type": "code", "execution_count": 12, "id": "e6fad82f", "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "ename": "LinAlgError", "evalue": "Singular matrix", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mLinAlgError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 2\u001b[0m \u001b[0mA\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m2.0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0mb\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1.0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m \u001b[0mx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinalg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msolve\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mA\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", "\u001b[0;32m<__array_function__ internals>\u001b[0m in \u001b[0;36msolve\u001b[0;34m(*args, **kwargs)\u001b[0m\n", "\u001b[0;32m~/my-conda-envs/jb/lib/python3.8/site-packages/numpy/linalg/linalg.py\u001b[0m in \u001b[0;36msolve\u001b[0;34m(a, b)\u001b[0m\n\u001b[1;32m 391\u001b[0m \u001b[0msignature\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m'DD->D'\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0misComplexType\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32melse\u001b[0m \u001b[0;34m'dd->d'\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 392\u001b[0m \u001b[0mextobj\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mget_linalg_error_extobj\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0m_raise_linalgerror_singular\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 393\u001b[0;31m \u001b[0mr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mgufunc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msignature\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0msignature\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mextobj\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mextobj\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 394\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 395\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mwrap\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mresult_t\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcopy\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mFalse\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m~/my-conda-envs/jb/lib/python3.8/site-packages/numpy/linalg/linalg.py\u001b[0m in \u001b[0;36m_raise_linalgerror_singular\u001b[0;34m(err, flag)\u001b[0m\n\u001b[1;32m 86\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 87\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m_raise_linalgerror_singular\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0merr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflag\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 88\u001b[0;31m \u001b[0;32mraise\u001b[0m \u001b[0mLinAlgError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Singular matrix\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 89\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 90\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m_raise_linalgerror_nonposdef\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0merr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflag\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mLinAlgError\u001b[0m: Singular matrix" ] } ], "source": [ "# Case with infinitely many solution\n", "A = np.array([[2.0, 2], [2, 2]])\n", "b = np.array([1.0, 1])\n", "x = np.linalg.solve(A, b)" ] }, { "cell_type": "code", "execution_count": 13, "id": "daa8b7a8", "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "ename": "LinAlgError", "evalue": "Singular matrix", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mLinAlgError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 2\u001b[0m \u001b[0mA\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m2.0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0mb\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1.0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m \u001b[0mx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinalg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msolve\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mA\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", "\u001b[0;32m<__array_function__ internals>\u001b[0m in \u001b[0;36msolve\u001b[0;34m(*args, **kwargs)\u001b[0m\n", "\u001b[0;32m~/my-conda-envs/jb/lib/python3.8/site-packages/numpy/linalg/linalg.py\u001b[0m in \u001b[0;36msolve\u001b[0;34m(a, b)\u001b[0m\n\u001b[1;32m 391\u001b[0m \u001b[0msignature\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m'DD->D'\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0misComplexType\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32melse\u001b[0m \u001b[0;34m'dd->d'\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 392\u001b[0m \u001b[0mextobj\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mget_linalg_error_extobj\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0m_raise_linalgerror_singular\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 393\u001b[0;31m \u001b[0mr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mgufunc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msignature\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0msignature\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mextobj\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mextobj\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 394\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 395\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mwrap\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mresult_t\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcopy\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mFalse\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m~/my-conda-envs/jb/lib/python3.8/site-packages/numpy/linalg/linalg.py\u001b[0m in \u001b[0;36m_raise_linalgerror_singular\u001b[0;34m(err, flag)\u001b[0m\n\u001b[1;32m 86\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 87\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m_raise_linalgerror_singular\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0merr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflag\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 88\u001b[0;31m \u001b[0;32mraise\u001b[0m \u001b[0mLinAlgError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Singular matrix\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 89\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 90\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m_raise_linalgerror_nonposdef\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0merr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflag\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mLinAlgError\u001b[0m: Singular matrix" ] } ], "source": [ "# Case without solution\n", "A = np.array([[2.0, 2], [2, 2]])\n", "b = np.array([1.0, 0])\n", "x = np.linalg.solve(A, b)" ] }, { "cell_type": "markdown", "id": "3f37eb72", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "A unique solution does not exist if and only if the [determinant](https://en.m.wikipedia.org/wiki/Determinant) of $\\mathbf{A}$ is $0$, in which case $\\mathbf{A}$ is called a singular matrix. For a $2$-by-$2$ matrix, the determinant is defined as follows:\n", "\n", "$$ \n", "\\begin{aligned}\n", "\\operatorname{det}(A) &:= \\left| \n", "\\begin{matrix}\n", "a_{00} & a_{01}\\\\\n", "a_{10} & a_{11}\n", "\\end{matrix}\n", "\\right|\\\\\n", "&= a_{00}\\times a_{11} - a_{01}\\times a_{10}.\n", "\\end{aligned}\n", "$$" ] }, { "cell_type": "markdown", "id": "b1bc6b34", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "For example, the first system has a unique solution because\n", "\n", "$$\n", "\\left|\n", "\\begin{matrix}\n", "2 & 2\\\\\n", "0 & 2\n", "\\end{matrix}\n", "\\right|\n", "= 2\\times 2 - 2\\times 0 = 4>0.\n", "$$\n", "The last two systems do not have unique solutions because\n", "\n", "$$\n", "\\left|\n", "\\begin{matrix}\n", "2 & 2\\\\\n", "2 & 2\n", "\\end{matrix}\n", "\\right|\n", "= 2\\times 2 - 2\\times 2 = 0.\n", "$$\n", "We can use the function `det` from `np.linalg` to compute the determinant as follows:" ] }, { "cell_type": "code", "execution_count": 14, "id": "1fe88134", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(4.0, 0.0)" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.linalg.det?\n", "np.linalg.det(np.array([[2.0, 2], [0, 2]])), np.linalg.det(np.array([[2.0, 2], [2, 2]]))" ] }, { "cell_type": "markdown", "id": "4959390d", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "**Exercise** Use the `det` and `solve` functions to assign `x` to the `numpy` array storing the solution of the following linear system if the solution is unique else `None`.\n", "\n", "$$\n", "\\begin{aligned}\n", "x_0 + 2 x_1 + 3x_2 &= 14\\\\\n", "2x_0 + x_1 + 2x_2 &= 10\\\\\n", "3 x_0 + 2x_1 + x_2 &= 10.\n", "\\end{aligned}\n", "$$" ] }, { "cell_type": "code", "execution_count": 15, "id": "8a4f7719", "metadata": { "deletable": false, "nbgrader": { "cell_type": "code", "checksum": "c08555513b5b41d8863530cc1d01888d", "grade": false, "grade_id": "linalg", "locked": false, "schema_version": 3, "solution": true, "task": false }, "slideshow": { "slide_type": "-" }, "tags": [ "remove-output" ] }, "outputs": [ { "ename": "NotImplementedError", "evalue": "", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mNotImplementedError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0;31m# YOUR CODE HERE\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0;32mraise\u001b[0m \u001b[0mNotImplementedError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 3\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mNotImplementedError\u001b[0m: " ] } ], "source": [ "# YOUR CODE HERE\n", "raise NotImplementedError()\n", "x" ] }, { "cell_type": "code", "execution_count": 16, "id": "57305acc", "metadata": { "code_folding": [], "deletable": false, "editable": false, "nbgrader": { "cell_type": "code", "checksum": "f560fce9950282f5284df26adf91772c", "grade": true, "grade_id": "test-linalg", "locked": true, "points": 1, "schema_version": 3, "solution": false, "task": false }, "slideshow": { "slide_type": "-" }, "tags": [ "hide-input", "remove-output" ] }, "outputs": [ { "ename": "AssertionError", "evalue": "", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mAssertionError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 2\u001b[0m \u001b[0;31m# As the main test must be hidden, you may want to verify your solution\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0;31m# as explained in the lecture using matrix multiplication.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m \u001b[0;32massert\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mndarray\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", "\u001b[0;31mAssertionError\u001b[0m: " ] } ], "source": [ "# tests\n", "# As the main test must be hidden, you may want to verify your solution\n", "# as explained in the lecture using matrix multiplication.\n", "assert isinstance(x, np.ndarray) and x.shape == (3,)" ] }, { "cell_type": "code", "execution_count": 17, "id": "851b36d6", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "code", "checksum": "824c6e9697bed8dce93e16c5fe9d852f", "grade": true, "grade_id": "h_test-linalg", "locked": true, "points": 1, "schema_version": 3, "solution": false, "task": false }, "tags": [ "remove-cell" ] }, "outputs": [], "source": [ "# hidden tests" ] }, { "cell_type": "markdown", "id": "2e8c1827", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "## Solving non-linear equations" ] }, { "cell_type": "markdown", "id": "95a1e1b5", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "Suppose we want to solve:\n", "\n", "$$\n", "f(x) = 0\n", "$$\n", "for some possibly non-linear real-valued function $f(x)$ in one real-valued variable $x$. Quadratic equation with an $x^2$ term is an example. The following is another example." ] }, { "cell_type": "code", "execution_count": 18, "id": "3636d660", "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEYCAYAAABfgk2GAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAAuLklEQVR4nO3deXxV1b338c+PTBASxkCYCQgODIoQRodqRYtcFeehdbaltrXt7b29HexzW++9nZ4O93FqtdTqcaBStChYkMkBREUIEKYAAiGQkARCgAQIIdN6/siBxpjhhHNyxu/79Tovzj57Za/fzoYvOzt7r2XOOUREJPp1CHUBIiISHAp8EZEYocAXEYkRCnwRkRihwBcRiREKfBGRGKHAFxGJEQp8EZEYocCXzzGzPDObGoR+zjOzDWZ2zMy+00ybXma2zMyOmNlfvJ/9ysz+1YftrzGzkQEuu6X+fKqrnWtocp/DuTYJHgV+jPKG+kkzO25mB8zsBTNLOYtt+PMfww+A951zqc65J5tp82Ngp3Ouu3PuITPrBdwL/MmH7f8O+G8/6vNZG+vyp59HzCzLzE6ZmaeJJp/b52DUZmZJZvYXM9vr/Q98g5ld21ptElwK/Nh2vXMuBRgLjAf+T5D7HwxsbaXNVOC1Bsv3A4uccyd92P4C4Eoz63t25bXJ/fhelz8KgZ8Dzzezvql9DkZt8UA+8AWgK/CfwFwzy2ilNgkiBb7gnNsPvA2MarzOzC4ws/fN7KiZbTWzG7yfvwwMAt7y/pTwg6a23cLXvwtcCTzt/fpzG31dopmVAaO9fWz2rroWWNGg3W/M7I0Gy781s3fMLME5VwmsA6452+9No5qa7atxXT60PyvOuXnOuTeB0mbWN7XPPn/P/KjrhHPuMedcnnOuzjn3D2APMK6V2iSI4kNdgISemQ0EpgPzGn2eALxF/dnkNcClwHwzy3TO3WNmlwFfdc4tb2a7LX39F83sfeAV59xzjb/WOVdlZpOB95xz6Q1WjQZ2NFj+v8BuMxsDTAKmAZc656q967cBFzVR2z+89TRllXPuuiY+b7YvM2tcly+1tZfG+9zW75nfzCwdOJfP/wTX5PGQ4FDgx7Y3zawGKAMWAr9stH4SkAL82jlXB7zrDcq7gMd82L6/Xz8G2Njos27AsdMLzrlSM3sceIn6SwmXOufKGrQ/BnzuEkIzgd6iVvr6TF0+1tZeGu/zZ2pr77q8/9HPBl50zm1vpTYJIl3SiW03Oue6OecGO+e+2cQ13n5AvjesT9sL9Pdx+/5+/Rg+H/hHgNRGn22g/iz2x865/EbrUoGjPvbni+b6aqqu1mrDe7nLNfNadZY1Nt7ntn7Pzro2M+sAvAxUAY/4UJsEkQJfWlIIDPT+Iz5tELDf+761yRRa+/rWXMTnA38T9ZcKAPBeSnkGeBF4sIltXNDENjCzt72/O2jq9XZTxbTS12fq8rE2nHNXOOesmVdzl5xa03if2/o9O6vazMyAvwDpwC3NXCJq8nhIcCjwpSWfACeAH5hZgpldAVwPzPGuPwAM9ePrW9NU4C+i/k4QzKw/9b8jeBj4JjDa2wfe9UnU/9JwWeMNO+eudc6lNPNqfDthq301rMvH9mfFzOLNrCMQB8SZWUczi2+wvql99vl75qdnqA/065u6I6il4yFB4pzTKwZfQB4wtbV1wEjq7/AoA3KAmxq0mwHso/5H9O83s62Wvv596n/p29TX9QFOAQmNPk8DCqi/9rwR+E6Ddd8HPmywfBswLwDfqy4+9HW6rk6+tPejlseo/8mq4euxlva5Ld8zP+oa7K2lEjje4PWVQB8Pvc7+Zd4DIRIxzOyXwEHn3OOttPsEeMg5tyWc6mrnGprc53CuTYJHgS8iEiN0DV9EJEYo8EVEYoQCX0QkRvj9pK33sfyXqL+rog6Y5Zx7olEbA56g/vH9CuB+59z61radlpbmMjIy/C1RRCRmrFu37pBzrldT6wIxtEIN8O/OufVmlgqsM7NlzrmcBm2uBYZ7XxOpv193YmsbzsjIICsrKwAliojEBjPb29w6vy/pOOeKTp+tO+eOUT84UuNH52cAL7l6q4FupiFSRUSCKqDX8K1+7OuLqX/CsqH+1I+VfVoBvo+nIiIiARCwwLf62ZL+Dvyrc6688eomvqTJBwDMbKbVz+iTVVJSEqjyRERiXkAC3zsc6t+B2c65eU00KQAGNlgeQP3AWp/jnJvlnMt0zmX26tXk7x1EROQs+B34DUbI2+ac+99mmi0A7rV6k4Ay51yRv32LiIjvAnGXziXAPcBmM8v2fvYo9cPg4px7lvrR+qYDu6i/LfOBAPQrIiJt4HfgO+dW0fQ1+oZtHPAtf/sSEZGzpydtRUTCyHvbD/L8qj1U19a13riNFPgiImHk+Q/38OLHecR3aPHCyVlR4IuIhImDxyr5cNchZlzUj/r7YQJLgS8iEiYWbiqizsENY/q1y/YV+CIiYWJ+diEj+3VhWO/Udtm+Al9EJAzsLT1Bdv5RZrTT2T0o8EVEwsKC7ELM4PqLFPgiIlHLOceb2fsZn9GDvl07tVs/CnwRkRDLKSpnd8mJdr2cAwp8EZGQW5BdSHwHY/qo9p0mRIEvIhJCdXWOBRsL+cK5vejeObFd+1Lgi4iE0Nq8wxSVVbbbvfcNKfBFREJo/sZCOiXEcfWI9HbvS4EvIhIiVTV1LNpcxDUj00lODMRo9S1T4IuIhMgHO0s4WlHd7nfnnKbAFxEJkfnZhXRPTuCy4cGZzlWBLyISAidO1bAs5wDTR/clIS44UazAFxEJgeXbDnCyupYZY/oHrc+ABL6ZPW9mB81sSzPrrzCzMjPL9r5+Goh+RUQi1Zsb9tOva0cyB3cPWp+BOsP3ANNaafOBc26M9/XfAepXRCTiFJdVsuLTEm68uD8d2mFmq+YEJPCdcyuBw4HYlohItPv7+gLqHNyeOTCo/QbzGv5kM9toZm+b2cjmGpnZTDPLMrOskpKSIJYnItL+6uocc7PymTS0BxlpnYPad7ACfz0w2Dl3EfAU8GZzDZ1zs5xzmc65zF69gnOrkohIsKzeU8re0gruGB/cs3sIUuA758qdc8e97xcBCWaWFoy+RUTCydy1+aR2jOfadh4ZsylBCXwz62PeKdjNbIK339Jg9C0iEi7KKqpZtKWYmy7uT8eEuKD3H5DBG8zsVeAKIM3MCoCfAQkAzrlngVuBb5hZDXASuNM55wLRt4hIpJi/cT9VNXVB/2XtaQEJfOfcXa2sfxp4OhB9iYhEqjlr8hnVvwuj+ncNSf960lZEJAi27C8jp6icO0J0dg8KfBGRoJizdh9J8R24IYhDKTSmwBcRaWcnq2qZv6GQ6aP70rVTQsjqUOCLiLSzt7cUcexUTUjuvW9IgS8i0s7mrM0no2cyE4f0CGkdCnwRkXaUW3KcNXsOc/v4gXgfRwoZBb6ISDuam1VAXAfj1rEDQl2KAl9EpL1U19bx9/UFXHleb3p36RjqchT4IiLtZdHmIkqOneIrEweFuhRAgS8i0m5e+DCPoWmd+cK54THyrwJfRKQdbNh3hOz8o9w3JSOos1q1RIEvItIOXvgwj9SkeG4ZF/pf1p6mwBcRCbDiskoWbS7i9vEDSUkKyBiVAaHAFxEJsFdW76XWOe6bnBHqUj5DgS8iEkCV1bX8dc0+rjo/nUE9k0Ndzmco8EVEAmjBxkIOn6jiwUsyQl3K5yjwRUQCxDnHCx/mcV56KpPP6Rnqcj4nIIFvZs+b2UEz29LMejOzJ81sl5ltMrOxgehXRCScfLLnMNuKynngkoyQj5vTlECd4XuAaS2svxYY7n3NBJ4JUL8iImHjhQ/30D05gRsvDt0kJy0JSOA751YCh1toMgN4ydVbDXQzs76B6FtEJBzkH65gWc4B7powiI4JcaEup0nBuobfH8hvsFzg/exzzGymmWWZWVZJSUlQihMR8ddLH+dhZtwzeXCoS2lWsAK/qYtZrqmGzrlZzrlM51xmr17hMf6EiEhLTpyqYc7afK4d1Ye+XTuFupxmBSvwC4CGc3sNAAqD1LeISLt6dc0+jlXW8MAlQ0JdSouCFfgLgHu9d+tMAsqcc0VB6ltEpN1UVtcya2Uuk4f2ZNzg7qEup0UBGeTBzF4FrgDSzKwA+BmQAOCcexZYBEwHdgEVwAOB6FdEJNRey8rn4LFTPH7nmFCX0qqABL5z7q5W1jvgW4HoS0QkXFTV1PHsilzGDe7O5KHh96BVY3rSVkTkLL2xoYD9R0/y7S8OC8sHrRpT4IuInIWa2jr++P5uRvfvGjYzWrVGgS8ichb+samIvaUVPBIhZ/egwBcRabO6OsfT7+3i/D6pXH1BeqjL8ZkCX0SkjRZvLWbXweN868phYTNfrS8U+CIibeCc46l3dzG0V2emj46sIcEU+CIibfDOtoNsKyrnW1cMIy6Czu5BgS8i4jPnHE+9t4uBPTpxw5h+oS6nzRT4IiI+WrnzEBvzj/LNK4aREBd58Rl5FYuIhEBdneM3i7fTv1snbh4bnhOctEaBLyLigzez97O1sJwfTDuPpPjwnOCkNQp8EZFWVFbX8rslOxjdvyvXXxh51+5PU+CLiLTihQ/zKCyr5NHpF0TUffeNKfBFRFpw+EQVf3xvF1ed35vJ54T/iJgtUeCLiLTgqXd3cqKqhh9de36oS/GbAl9EpBl7S0/wyuq93DF+IMPTU0Ndjt8U+CIizfjN4h0kxHXge1PPDXUpAaHAFxFpwvp9R1i4uYivXTaU3l06hrqcgAhI4JvZNDPbYWa7zOxHTay/wszKzCzb+/ppIPoVEWkPzjl+uXAbaSlJzLx8aKjLCRi/57Q1szjgD8DVQAGw1swWOOdyGjX9wDl3nb/9iYi0t6U5B8jae4Rf3DSKzkkBmfo7LATiDH8CsMs5l+ucqwLmADMCsF0RkaCrqKrhf/6Rw7DeKdyROTDU5QRUIAK/P5DfYLnA+1ljk81so5m9bWYjm9uYmc00sywzyyopKQlAeSIivnti+U4KjpzkFzeOIj4CB0hrSSD2pqnHzlyj5fXAYOfcRcBTwJvNbcw5N8s5l+mcy+zVKzImBhaR6LC1sIznVu3hjsyBTBwa2Q9ZNSUQgV8ANPy5ZwBQ2LCBc67cOXfc+34RkGBmaQHoW0QkIGrrHI/O20z35AR+PD3yH7JqSiACfy0w3MyGmFkicCewoGEDM+tj3mndzWyCt9/SAPQtIhIQL3+cx8aCMv7zuhF0S04MdTntwu9fPzvnaszsEWAJEAc875zbamYPe9c/C9wKfMPMaoCTwJ3OucaXfUREQqKo7CS/XbKDy4anccNFkTsaZmsCcr+R9zLNokafPdvg/dPA04HoS4Kvrs6x93AFOYXlbC8up+TYKcpOVlN2spryyvo/yyqqcUDnxHiSk+LonBhPp8Q4OifG0T05kcE9O5ORlszQtBQy0pJJ7ZgQ6t0SOeNn87dS6xy/uHE03osRUSl6bjCVgDlYXsn7O0rYUljG1sJytheVc6KqFoC4DkbPzol07ZRAl04J9E7tyPDeqXTpGI+ZUVFVQ0VVLRVVtZw4VcOh41XsKD7GvA37P9NHWkoiQ3ulcPGgbowf3INxg7vTvXN0/hgt4W3J1mKW5hzgh9POZ1DP5FCX064U+AJA/uEKlmwt5u0txazfdwTnICUpnhF9u3Bb5kBG9O3CiH5dGJ6eclaz/VRW17K3tII9h06QV3qCvEMn2HHgGM+v2sOfVuQCMKx3CpmDu5OZ0YPLz02jd2p0PM4u4etYZTU/m7+V8/uk8tXLhoS6nHanwI9hB8oreS0rn8Vbi9myvxyAEX278L2p5/KlkX0Y3jslYJM9dEyI47w+qZzX57MjDlZW17KpoIy1eYdZt/cIizYXMWdtPmYwZmA3rh6RzjUj0jmnV0pU/6gtofH7pZ9y4Fglz9w9NiInJW8rC+ffnWZmZrqsrKxQlxF1dpccZ9aKXN7YsJ+q2jouHtSNa0f14Usj+zC4Z+eQ1lZX59hefIzl2w6wLOcAm/eXATAkrTNXj0jnugv7Mrp/V4W/+G3FpyXc9/wa7p+SwWM3NPssaMQxs3XOucwm1ynwY8emgqM88/5uFm8tJjGuA7dnDuRrlw0N6+uWRWUnWZ5zgKU5B1idW0p1rWN47xRuHTeAmy7uHzWjGEpwlRw7xbVPrKRn5yTmP3IJHRMic1LypijwY9y6vYf532Wf8uGuUlI7xnPv5MHcP2UIvVKTQl1am5RVVPOPzYX8fV0B6/cdpYPBZcN7ceu4AVw9Ij2q/tFK+6mrc9z3whrW7DnMW9++lHOjYGKThloKfF3Dj2JlFdX8evF2Xl2zj96pSTw6/XzumjAoYm+J7JqcwFcmDuYrEwezu+Q489YXMG/9fr796ga6Jydwx/hB3DN5MP27dQp1qRLGnluVywc7D/GLm0ZFXdi3Rmf4Ucg5x4KNhfzPP3I4UlHNg5dk8L2rzyU5Mfr+f6+tc3y8u5RXVu9laU4xAFePSOe+KRlMHtpT1/rlMzbmH+WWZz5i6gXpPHP32Kj8+6Ez/Biyr7SC/zN/Cys/LeGiAV3xPDCBUf27hrqsdhPXwbh0eBqXDk9j/9GTvLJ6L3PW7GPJ1gOcl57KfVMyuHlsf13uEY5VVvOdORvonZrEr2+J7gesmqMz/CjhnOO5D/bwu6X1c3B+/5pzuWdyBnEBuq0yklRW17JgYyEvfpTH1sJy0lKSePDSDO6eNJguEXo5S/z3vb9lMz97P3/7+mTGZ/QIdTntRr+0jXLHT9XwH69t5O0txVw9Ip3/njGSvl11Hds5x8e5pTy7IpeVn5aQkhTPVyYN4qFLhujunhgzb30B/zZ3I9+bei7fnTo81OW0KwV+FMstOc7XX17H7pLjPDr9Ah66dEhM/qjami37y/jTylwWbiokvkMHbh7bn29eMSysb0mVwNhWVM4tz3zEqP5defVrk6L+p14FfpR6Z9sB/nVONgnxHXj6rouZMkxTDLRmb+kJ/vxBLnOzCqitc9x0cX8euXIYGWmhfeBM2sfB8kpu/MOH1DrH/G9dSp+u0f+TnQI/ytTVOZ58dyePL9/JqP5dePbucQzorjPVtjhYXsmzK3KZ/cleauocM8b045ErhzG0V0qoS5MAqaiq4Y4/rWZ3yXHmfn1yVN+80JACP4qcrKrlO3M2sCznADeP7c8vbxqtO1D8cPBYJbNW5PLKJ3upqqnjhov68Z2rhiv4I1xdneMbs9exLOcAs+7JZOqI9FCXFDQK/ChRUVXDQ54sVu8p5afXjeD+KRm6Xh8gJcdO8ecPcnn5472cqqnl5rED+O5VwxnYQz85RaJfLtrGrJW5/PS6ETx4afSPgtmQAj8KHD9Vw4MvrCVr72H+9/Yx3Hhx/1CXFJVKjp3i2RW7eXn1XurqHLePH8gjVw6jn57ejRh//WQfj76xmXsnD+a/bhgZcydFLQV+QMYDNbNpZrbDzHaZ2Y+aWG9m9qR3/SYzGxuIfmNFeWU19/7lE9btO8KTd12ssG9HvVKT+M/rRrDyP67krgmDeC0rnyt++z6PLdjKwWOVoS5PWvHBzhL+c/4WrjivFz+9bkTMhX1r/D7DN7M44FPgaqCA+knN73LO5TRoMx34NjAdmAg84Zyb2Nq2dYZfPx7Ovc9/Qk5ROU/dNZZpo/qEuqSYUnCkgqfe2cXr6wtIiDPum5LBw5efo9m5wtCW/WXcNWs1/bt34rWHJ0fsmFH+au8z/AnALudcrnOuCpgDzGjUZgbwkqu3GuhmZn1b23BpaSnZ2dkA1NbW4vF42LRpEwDV1dV4PB62bNkCQGVlJR6Ph23btgFQUVGBx+Nhx44dABw/fhyPx8OuXbsAKCsrw+PxkJtbP9vSkSNH8Hg85OXlAXDo0CE8Hg/5+fkAHDx4EI/Hw/799VP1FRcX4/F4KC6uH79l//79eDweDh48CEB+fj4ej4dDhw4BkJeXh8fj4ciRIwDk5ubi8XgoK6sf733Xrl14PB6OHz8OwI4dO3juLy9wz59Wsq3oGD//QneKsxZTWVl/lrllyxY8Hg/V1dUAbNq0CY/HQ21t/VSE2dnZeDyeM9/LdevW8dJLL51ZXrt2LbNnzz6zvHr1al599dUzyx999BFz5849s7xq1Spef/31M8srVqxg3rx5Z5bfe+895s+ff2Z5+fLlvPXWW2eWly5dysKFC88sL168mMWLF59ZXrhwIUuXLj2z/NZbb7F8+fIzy/Pnz+e99947szxv3jxWrFhxZvn1119n1apVZ5bnzp3LRx99dGb51VdfZfXq1WeWZ8+ezdq1a88sv/TSS6xbt+7MssfjITs7mwHdk/nlTSP5jyFF3DCgmlkrc7nyN8v5+f97hqwNG4Ho/Lvn8XioqKgAYNu2bXg8nrD+u/fCX+fy5T+vpkunBB4eXsH7y/75dytS/+7B2eVeSwIR+P2B/AbLBd7P2toGADObaWZZZpZ1+i9ULDp+qoaconL2lJ5g1r3juHBgt1CXFNM6JsRx9+QMFn/3ciYN7UnBkQoeW7CVP63YTaV3vl8JjaKySlZ+eoiuyQn87euTSEnSEGHNCcQlnduALznnvupdvgeY4Jz7doM2C4FfOedWeZffAX7gnFvX1DZPi9VLOqdqarn7uU/YWFDG8/eN59LheqAq3GwqOMrvln7Kyk9L6JWaxLeuOIe7Jg46q/l+5eytzTvM/c+voXeXjvz1axM1pAjtf0mnABjYYHkAUHgWbYT68V9++Pom1uYd4fe3XaSwD1MXDujGSw9OYO7XJzMkrTOPvZXDlb99n1fX7KO6ti7U5cWEj3eXcu9f1tCna0fmzJyksPdBIAJ/LTDczIaYWSJwJ7CgUZsFwL3eu3UmAWXOuaIA9B11nnhnJ29mF/L9a87l+ov6hbocacWEIT3428xJvPLQRNK7duTH8zZz1e9X8Pq6AmoU/O3mg50lPOBZw8AenZgzczLpGgzPJ34HvnOuBngEWAJsA+Y657aa2cNm9rC32SIgF9gF/Bn4pr/9RqM3N+zn8eU7uWXsAL515bBQlyM+Mqsfk3/eN6bw/P2ZpHaM5/uvbeSa/7eSNzYo+ANt8ZZiHnoxiyFpKbz6tUkRN1VnKOnBqzCxZs9h7n7uEy4e1I2XH5pIYnxAHpGQEHDOsWRrMY8v38n24mMMSevMI1cOY8aYfsTH6bierdo6x+PLP+Wpd3cxZmA3PA+Mp1uybo9tTE/ahrm8Qye46Y8f0j05kXnfnKK/xFGirs6xNOcAT76zk5yicjJ6JvOtK4dx08X9FfxtVFZRzXf/toH3d5RwR+ZA/mvGSI0h1QwFfhg7WlHFzX/8iCMVVbzxzUs0TG8Ucs6xLOcAT7yzk62F5QzqkczXvzCUW8YOUGj5YEfxMWa+nEXh0ZM8dsNIvjxhkJ6gbYECP0w553joxSxW7TzE7K9NjOpp16T+eL+z7SBPvbuTjQVlpKUk8cAl9VMvdu0Um0+FtmbhpiL+4/WNdE6K59m7xzJusP6NtEaTmIcpz0d5vLv9II9dP0JhHwPMjKkj0rnqgt5npl787ZIdPPP+br4ycRAPXjpEd5t4nThVw2+X7MDzUR7jBnfnma+M1bSUAaDAD5GthWX8atF2rjq/N/dNyQh1ORJEZsaUc9KYck7amakX//xBLi98mMf1F/Xj/ikZjB4QG5N1NGV5zgF+On8LhWWV3D8lg0enX6CbGAJEl3RCoKKqhuueWsWJUzW8/d3L6aGBuGLevtIKnluVy+vrCqioqmXsoG7cNyWDa0f1jZmwKy6r5L/e2srbW4o5Nz2FX908WpdwzoKu4YeZH7y+kdfWFTD7qxOZco6epJV/Kq+s5vWsAl76OI+80gp6pSbx5QmD+PLEQVF7uae2zvHK6r38dskOqmvr+O7U4Xz10qEx8x9doCnww8iCjYV859UNPHLlML7/pfNCXY6Eqbo6x4qdJbz4UR7v7yihg8Flw3txy7gBXDMiPSru7qmtcyzLKebp93axZX85lw1P4+c3jmJwT92p5g8FfpjIP1zB9Cc+YHh6CnO/Pln3YotP8g6d4PV1BcxbX0BhWSWpHeO57sJ+3DpuAGMHdYu4WxRP1dTy5ob9/GlFLrmHTpDRM5l/u+Y8rr+wb8TtSzhS4IeB6to6bnv2Y3aXHGfRdy7TXKnSZnV1jo9zS3l9XQFvbymisrqOQT2SuWZEOlePSCczowdxHcI3MI+fquHVT/bx3KpcDpSfYlT/LnzjC8OYNqpPWNcdaRT4YeB3S3bw9Hu7+MOXx/IvF7Y694tIi46fqmHR5iIWbiri492lVNXW0T05gS+eXx/+l5+bRnJi6G/CO1lVy4pPS1iytZjlOQc4dqqGKef05BtXnMOlw9J0Rt8OdB9+iG0rKueZFbu5ddwAhb0EREpSPLdnDuT2zIEcq6xm5aeHWJZTzLKcYv6+voDEuA6MHtCVzIzujB/cg3GDuwdtWsZjldW8u/0gi7cU8/6OEk5W19ItOYEvjerD3ZMGM0aT+YSMzvDbWW2d45ZnPiL/cAXv/PsXNE6OtKvq2jrW7jnM+5+WkJV3mM37y6iurf83Pqx3CuMGdWd4egpD0jqTkdaZgd2T/bob5lhlNduLj5FTWM7WwjJyisrZUXyM6lpH79QkvjSyD9NG9WHCkB4k6HdWQaEz/BCa/clesvOP8vgdYxT20u4S4jowZVgaU4bV3+5bWV3LpoIy1uYdZt3eIyzNKeZvWf+cOjSug9G/Wycy0jrTs3MiyYlxdE6Kp1NCHJ2T4khOjKe2zlF+spoy76u8sv7PorJK9pb+cw7VHp0TGdmvCw9dOpSpF/Rm7KDudNC1+bCiwG9HxWWV/GbxDi4bnsaMMZrMRIKvY0IcE4b0YMKQfz7AdOREFXtKT5B36AR7vK+9pRXsOXSck1W1nDhVy8nqz8/Tm5wYR5eOCXTtVP8a1a8rt40bwMh+XRnRrwu9U5N0TT7MKfDb0WMLtlJdW8fPbxylfwgSNrp3TqR750TGDurebJvaOsfJ6loqTtXQoYPRpWOCHoSKAgr8drIs5wCLtxbzg2nn6UESiThxHYyUpHhSkhQR0cSvo2lmPYC/ARlAHnC7c+5IE+3ygGNALVDT3C8UosXxUzX8dP4WzktP5WuXDQ11OSIigP9z2v4IeMc5Nxx4x7vcnCudc2OiPewBfr90B8Xllfzy5tG6M0FEwoa/aTQDeNH7/kXgRj+3F/E2FRzlxY/y+MrEQYwb3Pw1UhGRYPM38NOdc0UA3j97N9POAUvNbJ2ZzWxpg2Y208yyzCyrpKTEz/KCq67O8egbm0lLSeIH084PdTkiIp/R6jV8M1sO9Gli1U/a0M8lzrlCM+sNLDOz7c65lU01dM7NAmZB/YNXbegj5N7YsJ8t+8t54s4xdOmoKetEJLy0GvjOuanNrTOzA2bW1zlXZGZ9gYPNbKPQ++dBM3sDmAA0GfiRqrK6lt8v3cGFA7py/YW6515Ewo+/l3QWAPd5398HzG/cwMw6m1nq6ffANcAWP/sNO89/uIfCskoenX6Bni4UkbDkb+D/GrjazHYCV3uXMbN+ZrbI2yYdWGVmG4E1wELn3GI/+w0rh09U8cx7u5l6QW8mDe0Z6nJERJrk1334zrlS4KomPi8Epnvf5wIX+dNPuHvynZ2cqKrhh/pFrYiEMd0k7qe8Qyd4ZfVe7hg/iOHpqaEuR0SkWQp8P/1myXYS4zvwvauHh7oUEZEWKfD9sG7vERZtLmbm5UPpndox1OWIiLRIgX+WnHP8atE2eqUmabwcEYkICvyztGTrAbL2HuHfrj6XzhpRUEQigAL/LFTX1vGbxdsZ1juF28YNCHU5IiI+UeCfhfnZheQeOsEPp51PvEbDFJEIobRqo9o6xx/f28WIvl2YekFzY8WJiIQfBX4bLdxcRO6hE3z7i8M0baGIRBQFfhvU1Tn+8O4uhvdO4UsjmxpAVEQkfCnw22BpzgF2HDjGI18cpgHSRCTiKPB95Jzj6fd2ktEzmX8Z3TfU5YiItJkC30fv7yhhy/5yvnnlMN2ZIyIRScnlA+ccT767k/7dOnHTxf1DXY6IyFlR4Pvg492lbNh3lIevOIcEnd2LSIRSevngyXd3kt4lSU/VikhEU+C3Ym3eYVbnHmbm5efQMSEu1OWIiJw1BX4rnnp3Fz07J/LlCYNCXYqIiF/8Cnwzu83MtppZnZllttBumpntMLNdZvYjf/oMpk0FR1n5aQlfvWwonRJ1di8ikc3fM/wtwM3AyuYamFkc8AfgWmAEcJeZjfCz36D4y6o9pCTFc/cknd2LSOTzK/Cdc9uccztaaTYB2OWcy3XOVQFzgBn+9BsMB8orWbipiNsyB5DaMSHU5YiI+C0Y1/D7A/kNlgu8nzXJzGaaWZaZZZWUlLR7cc2ZvXovtc5x/5SMkNUgIhJIrU7VZGbLgaZGCvuJc26+D300NeiMa66xc24WMAsgMzOz2XbtqbK6ltmf7OOq83szuGfnUJQgIhJwrQa+c26qn30UAAMbLA8ACv3cZrt6a2MhpSeqeOCSIaEuRUQkYIJxSWctMNzMhphZInAnsCAI/Z4V5xwvfJjHeempTDmnZ6jLEREJGH9vy7zJzAqAycBCM1vi/byfmS0CcM7VAI8AS4BtwFzn3Fb/ym4/a/YcJqeonPsvydAEJyISVVq9pNMS59wbwBtNfF4ITG+wvAhY5E9fwfLCh3l0S07gxjEaJE1EoouetG0g/3AFS3OKuWvCID1oJSJRR4HfwMur92Jm3DNpcKhLEREJOAW+V0VVDXPW7GPayD7069Yp1OWIiAScAt9r3vr9lFfW8MAlGaEuRUSkXSjwqb8V0/NRHqP7d2Xc4O6hLkdEpF0o8IEPdh5i18HjPKBbMUUkiinwgVdW7yUtJZF/ubBvqEsREWk3MR/4B49V8u72g9wydgBJ8boVU0SiV8wH/rz1+6mpc9w+fmDrjUVEIlhMB75zjrlr8xmf0Z1zeqWEuhwRkXYV04G/Nu8IuYdOcMd4zWglItEvpgN/ztp9pCTFM310U8P9i4hEl5gN/PLKahZtLuKGMf1ITvRrDDkRkYgQs4G/ILuQyuo67tQva0UkRsRs4M/Nyuf8PqmM7t811KWIiARFTAZ+TmE5mwrKuHP8QD1ZKyIxIyYDf25WPonxHbjxYk1yIiKxw98pDm8zs61mVmdmmS20yzOzzWaWbWZZ/vTpr8rqWuatL2DayD50S04MZSkiIkHl7+0pW4CbgT/50PZK59whP/vz25KtxZRX1nCHflkrIjHG3zlttwERdR38b2vzGdijE5OH9gx1KSIiQRWsa/gOWGpm68xsZksNzWymmWWZWVZJSUlAi9hbeoKPdpdyR+ZAOnSInP+kREQCodUzfDNbDjT1KOpPnHPzfeznEudcoZn1BpaZ2Xbn3MqmGjrnZgGzADIzM52P2/fJa1kFdDC4dZwu54hI7Gk18J1zU/3txDlX6P3zoJm9AUwAmgz89lJX55i3voDLz+1Fn64dg9m1iEhYaPdLOmbW2cxST78HrqH+l71BlbX3CIVlldykWzFFJEb5e1vmTWZWAEwGFprZEu/n/cxskbdZOrDKzDYCa4CFzrnF/vR7NuZn76dTQhxTL0gPdtciImHB37t03gDeaOLzQmC6930ucJE//firuraORZuLuHpEOp2TNFCaiMSmmHjSdtXOQxypqOaGi/qFuhQRkZCJicCfn72frp0SuPzcXqEuRUQkZKI+8Cuqaliac4Dpo/uSGB/1uysi0qyoT8Dl2w5SUVXLjDG6nCMisS3qA39B9n76dOnIhIweoS5FRCSkojrwj5yo4v0dJdwwpp+GUhCRmBfVgf/2lmJq6pzuzhERIcoDf372fs7p1ZmR/bqEuhQRkZCL2sAvPHqSNXmHmTGmf0QN3ywi0l6iNvD/sakQ59DlHBERr6gN/PnZhVw0oCsZaZ1DXYqISFiIysDfdfAYWwvLuWGMRsYUETktKgN/QXYhZnD9hX1DXYqISNiIusB3zjF/YyFTzulJ7y6a6ERE5LSoGyv4ZHUtk4f2ZMqwtFCXIiISVqIu8JMT4/n1LReGugwRkbATdZd0RESkaf5OcfhbM9tuZpvM7A0z69ZMu2lmtsPMdpnZj/zpU0REzo6/Z/jLgFHOuQuBT4EfN25gZnHAH4BrgRHAXWY2ws9+RUSkjfwKfOfcUudcjXdxNTCgiWYTgF3OuVznXBUwB5jhT78iItJ2gbyG/yDwdhOf9wfyGywXeD8TEZEgavUuHTNbDvRpYtVPnHPzvW1+AtQAs5vaRBOfuRb6mwnMBBg0aFBr5YmIiI9aDXzn3NSW1pvZfcB1wFXOuaaCvAAY2GB5AFDYQn+zgFkAmZmZzf7HICIibePvXTrTgB8CNzjnKpppthYYbmZDzCwRuBNY4E+/IiLSdtb0SbmPX2y2C0gCSr0frXbOPWxm/YDnnHPTve2mA48DccDzzrlf+Lj9EmDvWZaXBhw6y68NN9GyL9GyH6B9CUfRsh/g374Mds71amqFX4EfzswsyzmXGeo6AiFa9iVa9gO0L+EoWvYD2m9f9KStiEiMUOCLiMSIaA78WaEuIICiZV+iZT9A+xKOomU/oJ32JWqv4YuIyGdF8xm+iIg0oMAXEYkRURP4ZtbDzJaZ2U7vn92baZdnZpvNLNvMsoJdZ3NaG0La6j3pXb/JzMaGok5f+LAvV5hZmfcYZJvZT0NRZ2vM7HkzO2hmW5pZH0nHpLV9iZRjMtDM3jOzbWa21cy+20SbiDguPu5LYI+Lcy4qXsBvgB953/8I+L/NtMsD0kJdb6Oa4oDdwFAgEdgIjGjUZjr1g9MZMAn4JNR1+7EvVwD/CHWtPuzL5cBYYEsz6yPimPi4L5FyTPoCY73vU6kflj1S/634si8BPS5Rc4ZP/ZDLL3rfvwjcGLpS2syXIaRnAC+5equBbmbWN9iF+iBqhsN2zq0EDrfQJFKOiS/7EhGcc0XOufXe98eAbXx+9N2IOC4+7ktARVPgpzvniqD+Gwn0bqadA5aa2TrvyJzhwJchpCNlmGlf65xsZhvN7G0zGxmc0gIuUo6JryLqmJhZBnAx8EmjVRF3XFrYFwjgcYmoScxbGqq5DZu5xDlXaGa9gWVmtt179hNKvgwh3aZhpkPIlzrXUz/ex3HvOEtvAsPbu7B2ECnHxBcRdUzMLAX4O/Cvzrnyxqub+JKwPS6t7EtAj0tEneE756Y650Y18ZoPHDj9Y5v3z4PNbKPQ++dB4A3qL0GEmi9DSLdpmOkQarVO51y5c+649/0iIMHM0oJXYsBEyjFpVSQdEzNLoD4gZzvn5jXRJGKOS2v7EujjElGB34oFwH3e9/cB8xs3MLPOZpZ6+j1wDdDkXQtB5ssQ0guAe713IEwCyk5fwgozre6LmfUxM/O+n0D938PSz20p/EXKMWlVpBwTb41/AbY55/63mWYRcVx82ZdAH5eIuqTTil8Dc83sIWAfcBuAfXao5nTgDe/3Lx74q3NucYjqPcM5V2NmjwBL+OcQ0lvN7GHv+meBRdTffbALqAAeCFW9LfFxX24FvmFmNcBJ4E7nvSUhnJjZq9TfJZFmZgXAz4AEiKxjAj7tS0QcE+AS4B5gs5llez97FBgEEXdcfNmXgB4XDa0gIhIjoumSjoiItECBLyISIxT4IiIxQoEvIhIjFPgiIjFCgS8iEiMU+CIiMUKBL9IG3vHLr/a+/7mZPRnqmkR8FU1P2ooEw8+A//YOvncxcEOI6xHxmZ60FWkjM1sBpABXeMcxF4kIuqQj0gZmNpr6mYpOKewl0ijwRXzkHXZ7NvUzKp0wsy+FuCSRNlHgi/jAzJKBecC/O+e2Af8DPBbSokTaSNfwRURihM7wRURihAJfRCRGKPBFRGKEAl9EJEYo8EVEYoQCX0QkRijwRURixP8HfcQVfXBnE0sAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "filenames": { "image/png": "/home/cs1302/www/cs1302book/_build/jupyter_execute/Lab9/Monte Carlo and Root Finding_40_0.png" }, "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "f = lambda x: x * (x - 1) * (x - 2)\n", "x = np.linspace(-0.5, 2.5)\n", "plt.plot(x, f(x))\n", "plt.axhline(color=\"gray\", linestyle=\":\")\n", "plt.xlabel(r\"$x$\")\n", "plt.title(r\"Plot of $f(x)=x(x-1)(x-2)$\")\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "f2c7eee0", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "While it is clear that the above function has three roots, namely, $x=0, 1, 2$, can we write a program to compute a root of any given continuous function $f$?" ] }, { "cell_type": "code", "execution_count": 19, "id": "1cbf921e", "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "data": { "text/html": [ "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%%html\n", "" ] }, { "cell_type": "markdown", "id": "191dd0cf", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "The following function `bisection` \n", "- takes as arguments \n", " - a continuous function `f`,\n", " - two real values `a` and `b`, \n", " - a positive integer `n` indicating the maximum depth of the recursion, and\n", "- returns a list `[xstart, xstop]` if the bisection succeeds in capturing a root in the interval `[xstart, xstop]` bounded by `a` and `b`, or else, returns a empty list `[]`." ] }, { "cell_type": "code", "execution_count": 20, "id": "8068038c", "metadata": { "code_folding": [ 13 ], "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "6ba8f476d2604502afc1d99d0b9de071", "version_major": 2, "version_minor": 0 }, "text/plain": [ "interactive(children=(FloatSlider(value=-0.5, description='a', max=2.5, min=-0.5, step=0.5), FloatSlider(value…" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEYCAYAAABfgk2GAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAAvn0lEQVR4nO3deXxV1bn/8c9DRjIQQkYIQ0AGZZIhjKJiBUWq4ohS59pS29rp3t4O9ndb79DhV9v+nFot1zpTkSoIKjI4AKIyBAxDwhwSCEkghJCEhJDknOf3xzncRszIOcmZnvfrlRdnWNnr2Wz4ZmedtdcWVcUYY0zw6+brAowxxnQNC3xjjAkRFvjGGBMiLPCNMSZEWOAbY0yIsMA3xpgQYYFvjDEhwgLfGGNChAW+8Usikisi04O1vwslIikiskZEKkTkbyLyWxH5YTu/d7OIjOjkEo0fE7vS1viCiBQAaYADaAA+BR5S1SNd1Pc3VPX9zu7L20TkT0C0qn5HRFKAHGCwqp5px/fOBe5Q1Vs7uUzjp+wM3/jSDaoaB/QGjgFP+bieQDAD+If78f3AivaEvdty4CoR6d0ZhRn/Z4FvfE5V64A3gOHnXhORAhGZ4X78UxE5KiLVIrJXRK52v95HRN4UkTIROSQi32+6XRHpJyJL3O+Xi8jTIvIK0B94W0ROi8hPmunvEhFZKyKn3EM9N5633QIR+bGI7BCRShF5XUSim9u31rbVwe1EikglMMpd+07gOmDdee1+LyJLmzx/TEQ+EJEI99/zVuCa1o6HCV4W+MbnRCQGuAPY2Mx7w4CHgQmqGg9cCxSISDfgbWA7kAFcDfxQRK51f18Y8A5QCGS62yxS1XuAw7h/u1DV35/XX4R7u6uBVOB7wEJ3HU3NBWYBA4HRuM62z6+9PdtqczsAqloPTAGOu+sehSv8957X9P/iOosfIyIPubd9i6o2uN/fDVzaXB8m+FngG196S0ROAVXATOCxZto4gChguPsstUBVDwITgBRV/U9VrVfVfOB/gDvd3zcR6AP8m6rWqGqdqm5oR02TgTjgd+7tfojrB8e889o9qarFqnoSV6iPucBttWc754zB9QPunJ5AddMGqloOPA68DPwcmK2qlU2aVLu/z4QgC3zjSzepak9cgf4wsE5E0ps2UNUDwA+BR4HjIrJIRPoAA4A+7qGSU+4fHI/g+iAYoB9QqKqNHaypD3BEVZ1NXivE9RtCU6VNHtfiCvYL2VZ7tnPOGL4Y+BVAfDPtPsd19v/zZj4EjwdOtdKHCWIW+MbnVNWhqktwnc1Pa+b9v6vqNFwhr7iGLY4Ah1S1Z5OveFWd7f62I0B/EQlvrstWyikG+rmHjM7pDxzt+J55dVvgGoppGvg7gKFNG4jIKOAZ4CXg681s45LztmFCiAW+8TlxmQMk4hpjbvreMBH5iohEAXXAGVw/GDYDVe4PdLuLSJiIjBSRCe5v3QyUAL8TkVgRiRaRy9zvHQMGtVDOJqAG+ImIRLjn5t8ALLqAXfPmtuDLgb8CuPLcExHJwDUs9BDwHWBU02sL3H+H44E1F9i/CXAW+MaX3haR07jG8H8N3Kequee1iQJ+B5zANfyRCjyiqg5c4TkGOOR+/zkgAVy/NbjfH4zrQ9oiXB8MA/wW+D/uoaAfN+3M/eHojbhmwJwA/gLcq6p7Orpz3tyWe6grEWj6vS8Ds90/8Hrg+gHwJ1Vdrqq1uD4T+XWT9jcCa1W1uKP9m+BgF14ZE8BE5De4Zu483o62m4AHVXVXpxdm/JIFvjHGhAgb0jHGmBBhgW+MMSHCAt8YY0JEc3OUO0RE+uGaLZAOOIEFqvrEeW0EeAKYjevikvtVdVtb205OTtbMzExPSzTGmJCxdevWE6qa0tx7Hgc+0Aj8q6puE5F4YKuIrFHVvCZtrgOGuL8m4bowZFJbG87MzCQ7O9sLJRpjTGgQkcKW3vN4SEdVS86dratqNa4LZ86/DH0O8LK6bAR62hKtxhjTtbw6hi8imcBYXFcYNpWB61L3c4r48g8FY4wxnchrgS8iccCbwA9Vter8t5v5lmYvABCR+SKSLSLZZWVl3irPGGNCnlcC373u95vAQvciWOcrwrV64Tl9cS0s9SWqukBVs1Q1KyWl2c8djDHGXACPA989A+dvwG5V/VMLzZYD97oXyZoMVKpqiad9G2OMaT9vzNK5DLgH2CkiOe7XHsG1DCyq+iyuRZ1mAwdwTct8wAv9GmOM6QCPA999F6HmxuibtlHgu572ZYwx5sLZlbbGGONHPtpznOc3HKLB4Wy7cQdZ4BtjjB95/pNDvPRZAeHdWh04uSAW+MYY4yeOV9fxyYETzLm0D675MN5lgW+MMX7i3R0lOBVuHNOnU7ZvgW+MMX5iWU4xI/r0YHBqfKds3wLfGGP8QGF5DTlHTjGnk87uwQLfGGP8wvKcYkTghkst8I0xJmipKm/lHGVCZi96J3TvtH4s8I0xxsfySqo4WFbTqcM5YIFvjDE+tzynmPBuwuyRnXubEAt8Y4zxIadTWb69mCuHppAYG9mpfVngG2OMD20pOElJZV2nzb1vygLfGGN8aNn2YrpHhDFzeFqn92WBb4wxPlLf6GTFzhKuGZFGTKQ3VqtvnQW+Mcb4yMf7yzhV29Dps3POscA3xhgfWZZTTGJMBJcP6ZrbuVrgG2OMD9ScbWRN3jFmj+pNRFjXRLEFvjHG+MD7u49xpsHBnDEZXdanVwJfRJ4XkeMisquF96eLSKWI5Li/fumNfo0xJlC99flR+iREkzUgscv69NYZ/ovArDbafKyqY9xf/+mlfo0xJuCUVtaxbl8ZN43NoFsn3NmqJV4JfFVdD5z0xraMMSbYvbmtCKfC3Kx+XdpvV47hTxGR7SLynoiMaKmRiMwXkWwRyS4rK+vC8owxpvM5ncri7CNMHtSLzOTYLu27qwJ/GzBAVS8FngLeaqmhqi5Q1SxVzUpJ6ZqpSsYY01U2HiqnsLyWOyZ07dk9dFHgq2qVqp52P14BRIhIclf0bYwx/mTxliPER4dzXSevjNmcLgl8EUkX9y3YRWSiu9/yrujbGGP8RWVtAyt2lXLz2AyiI8K6vH+vLN4gIq8B04FkESkCfgVEAKjqs8BtwLdFpBE4A9ypquqNvo0xJlAs236U+kZnl39Ye45XAl9V57Xx/tPA097oyxhjAtWizUcYmdGDkRkJPunfrrQ1xpgusOtoJXklVdzho7N7sMA3xpgusWjLYaLCu3FjFy6lcD4LfGOM6WRn6h0s+7yY2aN6k9A9wmd1WOAbY0wne29XCdVnG30y974pC3xjjOlki7YcITMphkkDe/m0Dgt8Y4zpRPllp9l86CRzJ/TDfTmSz1jgG2NMJ1qcXURYN+G2cX19XYoFvjHGdJYGh5M3txVx1bBUUntE+7ocC3xjjOksK3aWUFZ9lrsm9fd1KYAFvjHGdJoXPilgUHIsVw71j5V/LfCNMaYTfH64gpwjp7hvamaX3tWqNRb4xhjTCV74pID4qHBuHe/7D2vPscA3xhgvK62sY8XOEuZO6EdclFfWqPQKC3xjjPGyVzcW4lDlvimZvi7lCyzwjTHGi+oaHPx982GuvjiN/kkxvi7nCyzwjTHGi5ZvL+ZkTT1fvyzT16V8iQW+McZ4iarywicFDEuLZ8pFSb4u50u8Evgi8ryIHBeRXS28LyLypIgcEJEdIjLOG/0aY4w/2XToJLtLqnjgskyfr5vTHG+d4b8IzGrl/euAIe6v+cAzXurXGGP8xgufHCIxJoKbxvruJiet8Urgq+p64GQrTeYAL6vLRqCniPT2Rt/GGOMPjpysZU3eMeZN7E90RJivy2lWV43hZwBHmjwvcr/2JSIyX0SyRSS7rKysS4ozxhhPvfxZASLCPVMG+LqUFnVV4Dc3mKXNNVTVBaqapapZKSn+sf6EMca0puZsI4u2HOG6ken0Tuju63Ja1FWBXwQ0vbdXX6C4i/o2xphO9drmw1TXNfLAZQN9XUqruirwlwP3umfrTAYqVbWki/o2xphOU9fgYMH6fKYMSmL8gERfl9MqryzyICKvAdOBZBEpAn4FRACo6rPACmA2cACoBR7wRr/GGONr/8g+wvHqszx+5xhfl9ImrwS+qs5r430FvuuNvowxxl/UNzp5dl0+4wckMmWQ/11odT670tYYYy7Q0s+LOHrqDN/7ymC/vNDqfBb4xhhzARodTv6y9iCjMhL85o5WbbHAN8aYC/DOjhIKy2t5OEDO7sEC3xhjOszpVJ7+6AAXp8cz85I0X5fTbhb4xhjTQStzSzlw/DTfvWqw39yvtj0s8I0xpgNUlac+PMCglFhmjwqsJcEs8I0xpgM+2H2c3SVVfHf6YMIC6OweLPCNMabdVJWnPjpAv17duXFMH1+X02EW+MYY007r959g+5FTfGf6YCLCAi8+A69iY4zxAadT+f3KPWT07M4t4/zzBidtscA3xph2eCvnKLnFVfxk1jCiwv3zBidtscA3xpg21DU4+MOqvYzKSOCG0YE3dn+OBb4xxrThhU8KKK6s45HZlwTUvPvzWeAbY0wrTtbU85ePDnD1xalMucj/V8RsjQW+Mca04qkP91NT38jPrrvY16V4zALfGGNaUFhew6sbC7ljQj+GpMX7uhyPWeAbY0wLfr9yLxFh3fjRjKG+LsUrLPCNMaYZ2w5X8O7OEr55+SBSe0T7uhyv8Ergi8gsEdkrIgdE5GfNvD9dRCpFJMf99Utv9GuMMZ1BVfnNu7tJjoti/hWDfF2O13h8T1sRCQP+DMwEioAtIrJcVfPOa/qxql7vaX/GGNPZVucdI7uwgl/fPJLYKK/c+tsveOMMfyJwQFXzVbUeWATM8cJ2jTGmy9XWN/Jf7+QxODWOO7L6+bocr/JG4GcAR5o8L3K/dr4pIrJdRN4TkREtbUxE5otItohkl5WVeaE8Y4xpvyfe309RxRl+fdNIwgNwgbTWeGNvmrvsTM97vg0YoKqXAk8Bb7W0MVVdoKpZqpqVkhIYNwY2xgSH3OJKnttwiDuy+jFpUGBfZNUcbwR+EdD0956+QHHTBqpapaqn3Y9XABEikuyFvo0xxiscTuWRJTtJjIng57MD/yKr5ngj8LcAQ0RkoIhEAncCy5s2EJF0cd/WXUQmuvst90LfxhjjFa98VsD2okr+/frh9IyJ9HU5ncLjj59VtVFEHgZWAWHA86qaKyIPud9/FrgN+LaINAJngDtV9fxhH2OM8YmSyjM8tmovlw9J5sZLA3c1zLZ4Zb6Re5hmxXmvPdvk8dPA097oy3Q9p1MpPFlLXnEVe0qrKKs+S+WZBirPNFBV5/qzsrYBBWIjw4mJCiM2MpzukWHERoaRGBPJgKRYMpNjGJQcR2ZyDPHREb7eLWP+16+W5eJQ5dc3jcI9GBGUgmeCqfGa41V1rN1bxq7iSnKLq9hTUkVNvQOAsG5CUmwkCd0j6NE9gtT4aIakxtMjOhwRoba+kdp6B7X1DmrONnLidD17S6tZ8vnRL/SRHBfJoJQ4xvbvyYQBvRg/IJHE2OD8Ndr4t1W5pazOO8ZPZ11M/6QYX5fTqSzwDQBHTtayKreU93aVsu1wBaoQFxXO8N49uD2rH8N792B4nx4MSYu7oLv91DU4KCyv5dCJGgrKayg4UcPeY9U8v+EQf12XD8Dg1DiyBiSSldmLK4YmkxofHJezG/9VXdfAr5blcnF6PN+4fKCvy+l0Fvgh7FhVHf/IPsLK3FJ2Ha0CYHjvHvxoxlCuHZHOkNQ4r93sIToijGHp8QxL/+KKg3UNDnYUVbKl4CRbCytYsbOERVuOIAJj+vVk5vA0rhmexkUpcUH9q7bxjT+u3sex6jqeuXtcQN6UvKPEnz87zcrK0uzsbF+XEXQOlp1mwbp8ln5+lHqHk7H9e3LdyHSuHZHOgKRYn9bmdCp7Sqt5f/cx1uQdY+fRSgAGJscyc3ga14/uzaiMBAt/47F1+8q47/nN3D81k0dvbPFa0IAjIltVNavZ9yzwQ8eOolM8s/YgK3NLiQzrxtysfnzz8kF+PW5ZUnmG9/OOsTrvGBvzy2lwKENS47htfF9uHpsRNKsYmq5VVn2W655YT1JsFMsevozoiMC8KXlzLPBD3NbCk/xpzT4+OVBOfHQ4904ZwP1TB5ISH+Xr0jqksraBd3YW8+bWIrYdPkU3gcuHpHDb+L7MHJ4WVP9pTedxOpX7XtjM5kMneft70xgaBDc2aaq1wLcx/CBWWdvA71bu4bXNh0mNj+KR2Rczb2L/gJ0SmRATwV2TBnDXpAEcLDvNkm1FLNl2lO+99jmJMRHcMaE/90wZQEbP7r4u1fix5zbk8/H+E/z65pFBF/ZtsTP8IKSqLN9ezH+9k0dFbQNfvyyTH80cSkxk8P18dziVzw6W8+rGQlbnlQIwc3ga903NZMqgJBvrN1+w/cgpbn3mU2ZcksYzd48Lyn8fdoYfQg6X1/J/lu1i/b4yLu2bwIsPTGRkRoKvy+o0Yd2EaUOSmTYkmaOnzvDqxkIWbT7MqtxjDEuL576pmdwyLsOGewzVdQ18f9HnpMZH8btbg/sCq5bYGX6QUFWe+/gQf1jtugfnj68Zyj1TMgnz0rTKQFLX4GD59mJe+rSA3OIqkuOi+Pq0TO6ePIAeATqcZTz3o9dzWJZzlNe/NYUJmb18XU6nsQ9tg9zps4382z+2896uUmYOT+M/54ygd4KNY6sqn+WX8+y6fNbvKyMuKpy7JvfnwcsG2uyeELNkWxH/sng7P5oxlB/MGOLrcjqVBX4Qyy87zbde2crBstM8MvsSHpw2MCR/VW3LrqOV/HV9Pu/uKCa8WzduGZfBd6YP9uspqcY7dpdUcesznzIyI4HXvjk56H/rtcAPUh/sPsYPF+UQEd6Np+eNZepgu8VAWwrLa/ifj/NZnF2Ew6ncPDaDh68aTGayby84M53jeFUdN/35ExyqLPvuNNITgv83Owv8ION0Kk9+uJ/H39/PyIwePHv3ePom2plqRxyvquPZdfks3FRIo1OZM6YPD181mEEpcb4uzXhJbX0jd/x1IwfLTrP4W1OCevJCUxb4QeRMvYPvL/qcNXnHuGVcBr+5eZTNQPHA8eo6FqzL59VNhdQ3Ornx0j58/+ohFvwBzulUvr1wK2vyjrHgnixmDE/zdUldxgI/SNTWN/Lgi9lsPFTOL68fzv1TM2283kvKqs/yPx/n88pnhZxtdHDLuL784Ooh9OtlvzkFot+s2M2C9fn88vrhfH1a8K+C2ZQFfhA4fbaRr7+whezCk/xp7hhuGpvh65KCUln1WZ5dd5BXNhbidCpzJ/Tj4asG08eu3g0Yf990mEeW7uTeKQP4jxtHhNxJUWuB75X1QEVklojsFZEDIvKzZt4XEXnS/f4OERnnjX5b9eijnd5FV6mqa+Dev21i6+EKnpw3NnDDPgCOSUp8FP9+/XDW/9tVzJvYn39kH2H6Y2t5dHkux6vrfF2eacPH+8v492W7mD4shV9ePzzkwr4tHp/hi0gYsA+YCRThuqn5PFXNa9JmNvA9YDYwCXhCVSe1tW2PzvBFwI9/e2mvytoG7n1+E3klVTw1bxyzRqb7uqQLF4DHpKiilqc+OMAb24qICBPum5rJQ1dcZHfn8kO7jlYyb8FGMhK784+HpgTsmlGe6uylFSYCB1Q1393ZImAOkNekzRzgZfeNyzeKSE8R6a2qJa1tuLy8nJycHMaMGYPD4eCVV15h3LhxjB49moaGBhYuXEhWVhYjR46krq6ORYsWMWnSJC75618BKL34Ynr06EFMTAwOh4OysjISEhLo3r07jY2NnDhxgoSePekeHU1DYyPlJ07Qs2dPoqOjaWhooLy8nJ6JiURHRVHf0MDJ8nISExOJioqivr6ekydP0qtXLyIjIzl79iwVFRX0SkoiMiKCurNnOVVRQVJSEhEREdTV1XHq1CmSkpOJCA/nTF0dladOkZycTHh4OGfOnKGyspKUlBTCwsKora2lsrKKE44oHmlQBvYMx/lmLc7UVLp160ZNTQ3V1dWkpqXRTYTTNTWcrq4mLS0NEeH06dOcPn2a9HTXD4jq6mpqamtJT3N9eFVVXc2Z2lrSzj2vquJMXR1pqakAVFZVcfbsWVJTUlzPKyupr68nxf381KlTNDQ2kpLsmgpaceoUjsZGks89r6jA4XSSnJQEwMmKCppe27hy5UoAZs2aBcC7775LREQE11xzDQBvv/023bt3Z8aMGQAsW7aMHj16cNVVVwGwZMkSkpKSuPLKKwF44403SE9PZ9q0aQAsXryYvn37MnXqVABee+01Bg4cyOTJkwFYuHAhQ4cOZcKECQC8/PLLjBgxgvHjxwPw4osvMmbMGMaMGcNvbh7BwIot5DtTWLA+n9c3HuK2noXMmj6VrLGXfvHf3iWXUFtby+LFi5kyZQrDhg3j9OnTvPHGG0ybNo3BgwdTWVnJ0qVLueKKKxg0aBAVFRUsW7aM6dOnk5mZyYkTJ3jnnXe4+uqr6devH8ePH2fFihXMnDmTjIwMSktLWblyJbNmzSI9PZ2jR4+yZs0aZs+eTWpqKkeOHOGDDz7g+uuvJzk5mYKCAtauXcucOXNITEwkPz+f9evXc/PNN5OQkMCBAwfYsGEDt912G3Fxcezdu5fPPvuMuXPnEhMTw+7du9m0aRN33nkn0dHR7Nq1i+zsbO666y4iIiLYsWMH27Zt45577iEsLIycnBxycnK4//77Adi6dSu5ubnce++9AGzZsoV9+/Zx1113AbBx40YOHTrEvHnzAPj0008pKipi7ty5AGzYsIHS0lJuu+02ANatW0d5eTm33HILAB999BEFJWX8v30J9OgewUNDalm7ZiU33HADAKtXr6ahoYGvfvWrAfdvr0O55/631xpvDOlkAEeaPC9yv9bRNgCIyHwRyRaR7IaGBi+UF5gcqtTUN1LX4GRYWjyxUbbskS9FR4Rx95RMVv7gCiYPSqKoopZHl+fy13UHqXPf79f4RkllHev3nSAhJoLXvzWZOPu/0iJvDOncDlyrqt9wP78HmKiq32vS5l3gt6q6wf38A+Anqrq1tW2H6pDO2UYHdz+3ie1FlTx/3wSmDQmSC6oC+Jicb0fRKf6weh/r95WREh/Fd6dfxLxJ/S/ofr/mwm0pOMn9z28mtUc0f//mJFtShM7/0LYI6NfkeV+g+ALaGFzrv/z0jR1sKajgj7dfGjxhH2RG9+3Jy1+fyOJvTWFgciyPvp3HVY+t5bXNh2lwOH1dXkj47GA59/5tM+kJ0SyaP9nCvh28EfhbgCEiMlBEIoE7geXntVkO3OuerTMZqGxr/N5jv/pVp26+szzxwX7eyinmx9cM5YZL+/i6HO8K0GPSmokDe/H6/Mm8+uAk0hKi+fmSnVz9x3W8sbWIRgv+TvPx/jIeeHEz/Xp1Z9H8KaTZYnjt4pV5+O5ZOI8DYcDzqvprEXkIQFWfFdfcqKeBWUAt8ICqtjlWE2rz8N/6/Cg/fD2HW8f15Q+3j7YpZQFGVflo73H+uHofucVVDEqO5XtXD+aG0X0ID/PKDGgDrNxVyvcXfc5FKXG8+uBEkuIC61adnc0uvAoAmw+d5O7nNjG2f09eeXASkeEWEIFKVVmVW8rj7+9nT2k1A5NjefiqwcwZY8HvCYdTefz9fTz14QHG9OvJiw9MoGeMTY89nwW+nys4UcPNf/mExJhIlnxnqv0jDhJOp7I67xhPfrCfvJIqMpNi+O5Vg7l5bIYFfwdV1jbwg9c/Z+3eMu7I6sd/zBlha0i1wALfj52qreeWv3xKRW09S79zmS3TG4RUlTV5x3jig/3kFlfRv1cM37pyELeO62uh1Q57S6uZ/0o2xafO8OiNI/jaxP423NkKC3w/pao8+FI2G/afYOE3JwX1bdeM63h/sPs4T324n+1FlSTHRfHAZa5bLyZ0D82rQtvy7o4S/u2N7cRGhfPs3eMYP8D+j7TFbmLup178tIAP9xzn0RuGW9iHABFhxvA0rr4k9X9vvfjYqr08s/Ygd03qz9enDbTZJm41Zxt5bNVeXvy0gPEDEnnmrnF2W0ovsMD3kdziSn67Yg9XX5zKfVMzfV2O6UIiwtSLkpl6UfL/3nrxfz7O54VPCrjh0j7cPzWTUX1D42YdzXk/7xi/XLaL4so67p+aySOzL7FJDF5iQzo+UFvfyPVPbaDmbCPv/eAKetlCXCHvcHktz23I542tRdTWOxjXvyf3Tc3kupG9QybsSivr+I+3c3lvVylD0+L47S2jbAjnAtgYvp/5yRvb+cfWIhZ+YxJTL7Irac0/VdU18EZ2ES9/VkBBeS0p8VF8bWJ/vjapf9AO9zicyqsbC3ls1V4aHE5+MGMI35g2KGR+0HmbBb4fWb69mO+/9jkPXzWYH187zNflGD/ldCrr9pfx0qcFrN1bRjeBy4ekcOv4vlwzPC0oZvc4nMqavFKe/ugAu45WcfmQZP77ppEMSLKZap6wwPcTR07WMvuJjxmSFsfib02xudimXQpO1PDG1iKWbCuiuLKO+Ohwrh/dh9vG92Vc/54BN0XxbKODtz4/yl/X5ZN/oobMpBj+5Zph3DC6d8Dtiz+ywPcDDQ4ntz/7GQfLTrPi+5fbvVJNhzmdymf55byxtYj3dpVQ1+Ckf68YrhmexszhaWRl9iKsm/8G5umzjby26TDPbcjnWNVZRmb04NtXDmbWyHS/rjvQWOD7gT+s2svTHx3gz18bx1dH9/Z1OSbAnT7byIqdJby7o4TPDpZT73CSGBPBVy52hf8VQ5OJifT9JLwz9Q7W7StjVW4p7+cdo/psI1MvSuLb0y9i2uBkO6PvBDYP38d2l1TxzLqD3Da+r4W98Yq4qHDmZvVjblY/qusaWL/vBGvySlmTV8qb24qIDOvGqL4JZGUmMmFAL8YPSOyy2zJW1zXw4Z7jrNxVytq9ZZxpcNAzJoJrR6Zz9+QBjOnXs0vqMF9mZ/idzOFUbn3mU46crOWDf73S1skxnarB4WTLoZOs3VdGdsFJdh6tpMHh+j8+ODWO8f0TGZIWx8DkWDKTY+mXGOPRbJjqugb2lFaTV1xFbnEleSVV7C2tpsGhpMZHce2IdGaNTGfiwF5E2GdWXcLO8H1o4aZCco6c4vE7xljYm04XEdaNqYOTmTrYNd23rsHBjqJKthScZGthBavzSnk9+5+3Dg3rJmT07E5mcixJsZHERIYRGxVO94gwYqPCiIkMx+FUqs40UOn+qqpz/VlSWUdh+T/vodorNpIRfXrw4LRBzLgklXH9E+lmY/N+xQK/E5VW1vH7lXu5fEgyc8YE2c1MTECIjghj4sBeTBz4zwuYKmrqOVReQ8GJGg65vwrLazl04jRn6h3UnHVwpuHL9+mNiQyjR3QECd1dXyP7JHD7+L6M6JPA8D49SI2PsjF5P2eB34keXZ5Lg8PJf9800v4jGL+RGBtJYmwk4/onttjG4VTONDioPdtIt25Cj+gIuxAqCFjgd5I1ecdYmVvKT2YNswtJTMAJ6ybERYUTF2UREUw8Opoi0gt4HcgECoC5qlrRTLsCoBpwAI0tfaAQLE6fbeSXy3YxLC2eb14+yNflGGMM4PlNzH8GfKCqQ4AP3M9bcpWqjgn2sAf44+q9lFbV8ZtbRtnMBGOM3/A0jeYAL7kfvwTc5OH2At6OolO89GkBd03qz/gBLY+RGmNMV/M08NNUtQTA/WdqC+0UWC0iW0VkfmsbFJH5IpItItllZWUelte1nE7lkaU7SY6L4iezLvZ1OcYY8wVtjuGLyPtAejNv/aID/VymqsUikgqsEZE9qrq+uYaqugBYAK4LrzrQh88t/fwou45W8cSdY+gRbbesM8b4lzYDX1VntPSeiBwTkd6qWiIivYHjLWyj2P3ncRFZCkwEmg38QFXX4OCPq/cyum8CN4y2OffGGP/j6ZDOcuA+9+P7gGXnNxCRWBGJP/cYuAbY5WG/fuf5Tw5RXFnHI7MvsasLjTF+ydPA/x0wU0T2AzPdzxGRPiKywt0mDdggItuBzcC7qrrSw379ysmaep756CAzLkll8qAkX5djjDHN8mgevqqWA1c383oxMNv9OB+41JN+/N2TH+ynpr6Rn9oHtcYYP2aTxD1UcKKGVzcWcseE/gxJi/d1OcYY0yILfA/9ftUeIsO78aOZQ3xdijHGtMoC3wNbCytYsbOU+VcMIjU+2tflGGNMqyzwL5Cq8tsVu0mJj7L1cowxAcEC/wKtyj1GdmEF/zJzKLG2oqAxJgBY4F+ABoeT36/cw+DUOG4f39fX5RhjTLtY4F+AZTnF5J+o4aezLibcVsM0xgQIS6sOcjiVv3x0gOG9ezDjkpbWijPGGP9jgd9B7+4sIf9EDd/7ymC7baExJqBY4HeA06n8+cMDDEmN49oRzS0gaowx/ssCvwNW5x1j77FqHv7KYFsgzRgTcCzw20lVefqj/WQmxfDVUb19XY4xxnSYBX47rd1bxq6jVXznqsE2M8cYE5AsudpBVXnyw/1k9OzOzWMzfF2OMcZcEAv8dvjsYDmfHz7FQ9MvIsLO7o0xAcrSqx2e/HA/aT2i7KpaY0xAs8Bvw5aCk2zMP8n8Ky4iOiLM1+UYY8wFs8Bvw1MfHiApNpKvTezv61KMMcYjHgW+iNwuIrki4hSRrFbazRKRvSJyQER+5kmfXWlH0SnW7yvjG5cPonuknd0bYwKbp2f4u4BbgPUtNRCRMODPwHXAcGCeiAz3sN8u8bcNh4iLCufuyXZ2b4wJfB4FvqruVtW9bTSbCBxQ1XxVrQcWAXM86bcrHKuq490dJdye1Zf46Ahfl2OMMR7rijH8DOBIk+dF7teaJSLzRSRbRLLLyso6vbiWLNxYiEOV+6dm+qwGY4zxpjZv1SQi7wPNrRT2C1Vd1o4+mlt0RltqrKoLgAUAWVlZLbbrTHUNDhZuOszVF6cyICnWFyUYY4zXtRn4qjrDwz6KgH5NnvcFij3cZqd6e3sx5TX1PHDZQF+XYowxXtMVQzpbgCEiMlBEIoE7geVd0O8FUVVe+KSAYWnxTL0oydflGGOM13g6LfNmESkCpgDvisgq9+t9RGQFgKo2Ag8Dq4DdwGJVzfWs7M6z+dBJ8kqquP+yTLvBiTEmqLQ5pNMaVV0KLG3m9WJgdpPnK4AVnvTVVV74pICeMRHcNMYWSTPGBBe70raJIydrWZ1XyryJ/e1CK2NM0LHAb+KVjYWICPdMHuDrUowxxuss8N1q6xtZtPkws0ak06dnd1+XY4wxXmeB77Zk21Gq6hp54LJMX5dijDGdwgIf11TMFz8tYFRGAuMHJPq6HGOM6RQW+MDH+09w4PhpHrCpmMaYIGaBD7y6sZDkuEi+Orq3r0sxxphOE/KBf7y6jg/3HOfWcX2JCrepmMaY4BXygb9k21EancrcCf3abmyMMQEspANfVVm85QgTMhO5KCXO1+UYY0ynCunA31JQQf6JGu6YYHe0MsYEv5AO/EVbDhMXFc7sUc0t92+MMcElZAO/qq6BFTtLuHFMH2IiPVpDzhhjAkLIBv7ynGLqGpzcaR/WGmNCRMgG/uLsI1ycHs+ojARfl2KMMV0iJAM/r7iKHUWV3Dmhn11Za4wJGSEZ+IuzjxAZ3o2bxtpNTowxocPTWxzeLiK5IuIUkaxW2hWIyE4RyRGRbE/69FRdg4Ml24qYNSKdnjGRvizFGGO6lKfTU3YBtwB/bUfbq1T1hIf9eWxVbilVdY3cYR/WGmNCjKf3tN0NBNQ4+OtbjtCvV3emDErydSnGGNOlumoMX4HVIrJVROa31lBE5otItohkl5WVebWIwvIaPj1Yzh1Z/ejWLXB+SBljjDe0eYYvIu8DzV2K+gtVXdbOfi5T1WIRSQXWiMgeVV3fXENVXQAsAMjKytJ2br9d/pFdRDeB28bbcI4xJvS0GfiqOsPTTlS12P3ncRFZCkwEmg38zuJ0Kku2FXHF0BTSE6K7smtjjPELnT6kIyKxIhJ/7jFwDa4Pe7tUdmEFxZV13GxTMY0xIcrTaZk3i0gRMAV4V0RWuV/vIyIr3M3SgA0ish3YDLyrqis96fdCLMs5SveIMGZcktbVXRtjjF/wdJbOUmBpM68XA7Pdj/OBSz3px1MNDicrdpYwc3gasVG2UJoxJjSFxJW2G/afoKK2gRsv7ePrUowxxmdCIvCX5RwloXsEVwxN8XUpxhjjM0Ef+LX1jazOO8bsUb2JDA/63TXGmBYFfQK+v/s4tfUO5oyx4RxjTGgL+sBfnnOU9B7RTMzs5etSjDHGp4I68Ctq6lm7t4wbx/SxpRSMMSEvqAP/vV2lNDrVZucYYwxBHvjLco5yUUosI/r08HUpxhjjc0Eb+MWnzrC54CRzxmQE1PLNxhjTWYI28N/ZUYwqNpxjjDFuQRv4y3KKubRvApnJsb4uxRhj/EJQBv6B49XkFldx4xhbGdMYY84JysBfnlOMCNwwurevSzHGGL8RdIGvqizbXszUi5JI7WE3OjHGmHOCbq3gMw0OpgxKYurgZF+XYowxfiXoAj8mMpzf3Tra12UYY4zfCbohHWOMMc3z9BaHj4nIHhHZISJLRaRnC+1micheETkgIj/zpE9jjDEXxtMz/DXASFUdDewDfn5+AxEJA/4MXAcMB+aJyHAP+zXGGNNBHgW+qq5W1Ub3041A32aaTQQOqGq+qtYDi4A5nvRrjDGm47w5hv914L1mXs8AjjR5XuR+zRhjTBdqc5aOiLwPpDfz1i9UdZm7zS+ARmBhc5to5jVtpb/5wHyA/v37t1WeMcaYdmoz8FV1Rmvvi8h9wPXA1araXJAXAf2aPO8LFLfS3wJgAUBWVlaLPxiMMcZ0jKezdGYBPwVuVNXaFpptAYaIyEARiQTuBJZ70q8xxpiOk+ZPytv5zSIHgCig3P3SRlV9SET6AM+p6mx3u9nA40AY8Lyq/rqd2y8DCi+wvGTgxAV+r78Jln0Jlv0A2xd/FCz7AZ7tywBVTWnuDY8C35+JSLaqZvm6Dm8Iln0Jlv0A2xd/FCz7AZ23L3alrTHGhAgLfGOMCRHBHPgLfF2AFwXLvgTLfoDtiz8Klv2ATtqXoB3DN8YY80XBfIZvjDGmCQt8Y4wJEUET+CLSS0TWiMh+95+JLbQrEJGdIpIjItldXWdL2lpCWlyedL+/Q0TG+aLO9mjHvkwXkUr3McgRkV/6os62iMjzInJcRHa18H4gHZO29iVQjkk/EflIRHaLSK6I/KCZNgFxXNq5L949LqoaFF/A74GfuR//DPi/LbQrAJJ9Xe95NYUBB4FBQCSwHRh+XpvZuBanE2AysMnXdXuwL9OBd3xdazv25QpgHLCrhfcD4pi0c18C5Zj0Bsa5H8fjWpY9UP+vtGdfvHpcguYMH9eSyy+5H78E3OS7UjqsPUtIzwFeVpeNQE8R6d3VhbZD0CyHrarrgZOtNAmUY9KefQkIqlqiqtvcj6uB3Xx59d2AOC7t3BevCqbAT1PVEnD9RQKpLbRTYLWIbHWvzOkP2rOEdKAsM93eOqeIyHYReU9ERnRNaV4XKMekvQLqmIhIJjAW2HTeWwF3XFrZF/DicQmom5i3tlRzBzZzmaoWi0gqsEZE9rjPfnypPUtId2iZaR9qT53bcK33cdq9ztJbwJDOLqwTBMoxaY+AOiYiEge8CfxQVavOf7uZb/Hb49LGvnj1uATUGb6qzlDVkc18LQOOnfu1zf3n8Ra2Uez+8ziwFNcQhK+1ZwnpDi0z7UNt1qmqVap62v14BRAhIsldV6LXBMoxaVMgHRMRicAVkAtVdUkzTQLmuLS1L94+LgEV+G1YDtznfnwfsOz8BiISKyLx5x4D1wDNzlroYu1ZQno5cK97BsJkoPLcEJafaXNfRCRdRMT9eCKuf4flX9qS/wuUY9KmQDkm7hr/BuxW1T+10Cwgjkt79sXbxyWghnTa8DtgsYg8CBwGbgeQLy7VnAYsdf/9hQN/V9WVPqr3f6lqo4g8DKzin0tI54rIQ+73nwVW4Jp9cACoBR7wVb2taee+3AZ8W0QagTPAneqekuBPROQ1XLMkkkWkCPgVEAGBdUygXfsSEMcEuAy4B9gpIjnu1x4B+kPAHZf27ItXj4strWCMMSEimIZ0jDHGtMIC3xhjQoQFvjHGhAgLfGOMCREW+MYYEyIs8I0xJkRY4BtjTIiwwDemA9zrl890P/5vEXnS1zUZ017BdKWtMV3hV8B/uhffGwvc6ON6jGk3u9LWmA4SkXVAHDDdvY65MQHBhnSM6QARGYXrTkVnLexNoLHAN6ad3MtuL8R1R6UaEbnWxyUZ0yEW+Ma0g4jEAEuAf1XV3cB/AY/6tChjOsjG8I0xJkTYGb4xxoQIC3xjjAkRFvjGGBMiLPCNMSZEWOAbY0yIsMA3xpgQYYFvjDEh4v8DYmbAYFFqo5oAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "filenames": { "image/png": "/home/cs1302/www/cs1302book/_build/jupyter_execute/Lab9/Monte Carlo and Root Finding_44_1.png" }, "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "def bisection(f, a, b, n=10):\n", " if f(a) * f(b) > 0:\n", " return [] # because f(x) may not have a root between x=a and x=b\n", " elif n <= 0: # base case when recursion cannot go any deeper\n", " return [a, b] if a <= b else [b, a]\n", " else:\n", " c = (a + b) / 2 # bisect the interval between a and b\n", " return bisection(f, a, c, n - 1) or bisection(f, c, b, n - 1) # recursion\n", "\n", "\n", "# bisection solver\n", "import ipywidgets as widgets\n", "\n", "\n", "@widgets.interact(a=(-0.5, 2.5, 0.5), b=(-0.5, 2.5, 0.5), n=(0, 10, 1))\n", "def bisection_solver(a=-0.5, b=0.5, n=0):\n", " x = np.linspace(-0.5, 2.5)\n", " plt.plot(x, f(x))\n", " plt.axhline(color=\"gray\", linestyle=\":\")\n", " plt.xlabel(r\"$x$\")\n", " plt.title(r\"Bisection on $f(x)$\")\n", " [xstart, xstop] = bisection(f, a, b, n)\n", " plt.plot([xstart, xstop], [0, 0], \"r|-\")\n", " print(\"Interval: \", [xstart, xstop])" ] }, { "cell_type": "markdown", "id": "2311f74d", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "Try setting the values of $a$ and $b$ as follows and change $n$ to see the change of the interval step-by-step." ] }, { "cell_type": "code", "execution_count": 21, "id": "701c974f", "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/plain": [ "([-0.0009765625, 0.0], [1.0, 1.0009765625], [1.9998046875000002, 2.00234375])" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "bisection(f, -0.5, 0.5), bisection(f, 1.5, 0.5), bisection(f, -0.1, 2.5)" ] }, { "cell_type": "markdown", "id": "971fd477", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "**Exercise** Modify the function `bisection` to \n", "- take the floating point parameter `tol` instead of `n`, and\n", "- return the interval from the bisection method represented by a list `[xstart,xstop]` but as soon as the gap `xstop - xstart` is $\\leq$ `tol`." ] }, { "cell_type": "code", "execution_count": 22, "id": "be322724", "metadata": { "deletable": false, "nbgrader": { "cell_type": "code", "checksum": "6a4aae6b76256c6f757465fd4a93bed5", "grade": false, "grade_id": "bisection", "locked": false, "schema_version": 3, "solution": true, "task": false }, "slideshow": { "slide_type": "-" }, "tags": [ "remove-output" ] }, "outputs": [], "source": [ "def bisection(f, a, b, tol=1e-9):\n", " # YOUR CODE HERE\n", " raise NotImplementedError()" ] }, { "cell_type": "code", "execution_count": 23, "id": "8fd05bcd", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "code", "checksum": "a47780044fab40faea56393c61f7dcd9", "grade": true, "grade_id": "test-bisection", "locked": true, "points": 1, "schema_version": 3, "solution": false, "task": false }, "slideshow": { "slide_type": "-" }, "tags": [ "hide-input", "remove-output" ] }, "outputs": [ { "ename": "NotImplementedError", "evalue": "", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mNotImplementedError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 4\u001b[0m \u001b[0mf\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mlambda\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mx\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 5\u001b[0;31m \u001b[0mbisection\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1.5\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0.5\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 6\u001b[0m \u001b[0;32massert\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0misclose\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mbisection\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m-\u001b[0m\u001b[0;36m0.5\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0.5\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m9.313225746154785e-10\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0.0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mall\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 7\u001b[0m \u001b[0m_\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mbisection\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1.5\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0.5\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1e-2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m\u001b[0m in \u001b[0;36mbisection\u001b[0;34m(f, a, b, tol)\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mbisection\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtol\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1e-9\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2\u001b[0m \u001b[0;31m# YOUR CODE HERE\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 3\u001b[0;31m \u001b[0;32mraise\u001b[0m \u001b[0mNotImplementedError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", "\u001b[0;31mNotImplementedError\u001b[0m: " ] } ], "source": [ "# tests\n", "import numpy as np\n", "\n", "f = lambda x: x * (x - 1) * (x - 2)\n", "bisection(f, 1.5, 0.5)\n", "assert np.isclose(bisection(f, -0.5, 0.5), [-9.313225746154785e-10, 0.0]).all()\n", "_ = bisection(f, 1.5, 0.5, 1e-2)\n", "assert np.isclose(_, [1.0, 1.0078125]).all() or np.isclose(_, [0.9921875, 1.0]).all()\n", "assert np.isclose(\n", " bisection(f, -0.1, 2.5, 1e-3), [1.9998046875000002, 2.0004394531250003]\n", ").all()" ] }, { "cell_type": "code", "execution_count": 24, "id": "9cdeddf4", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "code", "checksum": "8160b43badb41f359045b5d097afe543", "grade": true, "grade_id": "h_test-bisection", "locked": true, "points": 1, "schema_version": 3, "solution": false, "task": false }, "tags": [ "remove-cell" ] }, "outputs": [], "source": [ "# hidden tests" ] } ], "metadata": { "jupytext": { "text_representation": { "extension": ".md", "format_name": "myst", "format_version": 0.13, "jupytext_version": "1.10.3" } }, "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.8" }, "source_map": [ 14, 18, 23, 27, 37, 41, 45, 49, 118, 128, 132, 137, 167, 172, 176, 192, 196, 200, 215, 220, 230, 236, 260, 285, 304, 308, 320, 331, 353, 357, 368, 379, 395, 421, 426, 438, 459, 484, 503, 507, 516, 534, 538, 547, 556, 588, 592, 600, 606, 627, 658 ], "widgets": { "application/vnd.jupyter.widget-state+json": { "state": { "153da0a246fc418f88cfaa4e01427558": { "model_module": "@jupyter-widgets/controls", "model_module_version": "1.5.0", "model_name": "FloatSliderModel", "state": { "_dom_classes": [], "_model_module": "@jupyter-widgets/controls", "_model_module_version": "1.5.0", "_model_name": "FloatSliderModel", "_view_count": null, "_view_module": "@jupyter-widgets/controls", "_view_module_version": "1.5.0", "_view_name": "FloatSliderView", "continuous_update": true, "description": "b", "description_tooltip": null, "disabled": false, "layout": "IPY_MODEL_a31d9848ee2a466fa3b8ee92ea83e9d2", "max": 2.5, "min": -0.5, "orientation": "horizontal", "readout": true, "readout_format": ".2f", "step": 0.5, "style": "IPY_MODEL_2482b59724fb4409b2dbab39be124485", "value": 0.5 } }, "2482b59724fb4409b2dbab39be124485": { "model_module": "@jupyter-widgets/controls", "model_module_version": "1.5.0", "model_name": "SliderStyleModel", "state": { "_model_module": "@jupyter-widgets/controls", "_model_module_version": "1.5.0", "_model_name": "SliderStyleModel", "_view_count": null, "_view_module": "@jupyter-widgets/base", "_view_module_version": "1.2.0", "_view_name": "StyleView", "description_width": "", "handle_color": null } }, "31689912f00e4837ab005ec2507a9a5d": { "model_module": "@jupyter-widgets/output", "model_module_version": "1.0.0", "model_name": "OutputModel", "state": { "_dom_classes": [], "_model_module": "@jupyter-widgets/output", "_model_module_version": "1.0.0", "_model_name": "OutputModel", "_view_count": null, "_view_module": "@jupyter-widgets/output", "_view_module_version": "1.0.0", "_view_name": "OutputView", "layout": "IPY_MODEL_40cbaa58a2c5461e8eab242344328294", "msg_id": "", "outputs": [ { "name": "stdout", "output_type": "stream", "text": "Interval: [-0.5, 0.5]\n" } ] } }, "3382016e18834c45b5b14ea3e427bdf1": { "model_module": "@jupyter-widgets/base", "model_module_version": "1.2.0", "model_name": "LayoutModel", "state": { "_model_module": "@jupyter-widgets/base", "_model_module_version": "1.2.0", "_model_name": "LayoutModel", "_view_count": null, "_view_module": "@jupyter-widgets/base", "_view_module_version": "1.2.0", "_view_name": "LayoutView", "align_content": null, "align_items": null, "align_self": null, "border": null, "bottom": null, "display": null, "flex": null, "flex_flow": null, "grid_area": null, "grid_auto_columns": null, "grid_auto_flow": null, "grid_auto_rows": null, "grid_column": null, "grid_gap": null, "grid_row": null, "grid_template_areas": null, "grid_template_columns": null, "grid_template_rows": null, "height": null, "justify_content": null, "justify_items": null, "left": null, "margin": null, "max_height": null, "max_width": null, "min_height": null, "min_width": null, "object_fit": null, "object_position": null, "order": null, "overflow": null, "overflow_x": null, "overflow_y": null, "padding": null, "right": null, "top": null, "visibility": null, "width": null } }, "40cbaa58a2c5461e8eab242344328294": { "model_module": "@jupyter-widgets/base", "model_module_version": "1.2.0", "model_name": "LayoutModel", "state": { "_model_module": "@jupyter-widgets/base", "_model_module_version": "1.2.0", "_model_name": "LayoutModel", "_view_count": null, "_view_module": "@jupyter-widgets/base", "_view_module_version": "1.2.0", "_view_name": "LayoutView", "align_content": null, "align_items": null, "align_self": null, "border": null, "bottom": null, "display": null, "flex": null, "flex_flow": null, "grid_area": null, "grid_auto_columns": null, "grid_auto_flow": null, "grid_auto_rows": null, "grid_column": null, "grid_gap": null, "grid_row": null, "grid_template_areas": null, "grid_template_columns": null, "grid_template_rows": null, "height": null, "justify_content": null, "justify_items": null, "left": null, "margin": null, "max_height": null, "max_width": null, "min_height": null, "min_width": null, "object_fit": null, "object_position": null, "order": null, "overflow": null, "overflow_x": null, "overflow_y": null, "padding": null, "right": null, "top": null, "visibility": null, "width": null } }, "450728e438624622b656eb74f2cdc46c": { "model_module": "@jupyter-widgets/controls", "model_module_version": "1.5.0", "model_name": "FloatSliderModel", "state": { "_dom_classes": [], "_model_module": "@jupyter-widgets/controls", "_model_module_version": "1.5.0", "_model_name": "FloatSliderModel", "_view_count": null, "_view_module": "@jupyter-widgets/controls", "_view_module_version": "1.5.0", "_view_name": "FloatSliderView", "continuous_update": true, "description": "a", "description_tooltip": null, "disabled": false, "layout": "IPY_MODEL_82a1713b423544e99b5bdf756111f0df", "max": 2.5, "min": -0.5, "orientation": "horizontal", "readout": true, "readout_format": ".2f", "step": 0.5, "style": "IPY_MODEL_5a445c3c64364598920c482ba07ace45", "value": -0.5 } }, "5a445c3c64364598920c482ba07ace45": { "model_module": "@jupyter-widgets/controls", "model_module_version": "1.5.0", "model_name": "SliderStyleModel", "state": { "_model_module": "@jupyter-widgets/controls", "_model_module_version": "1.5.0", "_model_name": "SliderStyleModel", "_view_count": null, "_view_module": "@jupyter-widgets/base", "_view_module_version": "1.2.0", "_view_name": "StyleView", "description_width": "", "handle_color": null } }, "6ba8f476d2604502afc1d99d0b9de071": { "model_module": "@jupyter-widgets/controls", "model_module_version": "1.5.0", "model_name": "VBoxModel", "state": { "_dom_classes": [ "widget-interact" ], "_model_module": "@jupyter-widgets/controls", "_model_module_version": "1.5.0", "_model_name": "VBoxModel", "_view_count": null, "_view_module": "@jupyter-widgets/controls", "_view_module_version": "1.5.0", "_view_name": "VBoxView", "box_style": "", "children": [ "IPY_MODEL_450728e438624622b656eb74f2cdc46c", "IPY_MODEL_153da0a246fc418f88cfaa4e01427558", "IPY_MODEL_ed302fda8d554243a2ffa1c0decee4af", "IPY_MODEL_31689912f00e4837ab005ec2507a9a5d" ], "layout": "IPY_MODEL_e2d07d875ad34e97a02d11bcacd7c417" } }, "82a1713b423544e99b5bdf756111f0df": { "model_module": "@jupyter-widgets/base", "model_module_version": "1.2.0", "model_name": "LayoutModel", "state": { "_model_module": "@jupyter-widgets/base", "_model_module_version": "1.2.0", "_model_name": "LayoutModel", "_view_count": null, "_view_module": "@jupyter-widgets/base", "_view_module_version": "1.2.0", "_view_name": "LayoutView", "align_content": null, "align_items": null, "align_self": null, "border": null, "bottom": null, "display": null, "flex": null, "flex_flow": null, "grid_area": null, "grid_auto_columns": null, "grid_auto_flow": null, "grid_auto_rows": null, "grid_column": null, "grid_gap": null, "grid_row": null, "grid_template_areas": null, "grid_template_columns": null, "grid_template_rows": null, "height": null, "justify_content": null, "justify_items": null, "left": null, "margin": null, "max_height": null, "max_width": null, "min_height": null, "min_width": null, "object_fit": null, "object_position": null, "order": null, "overflow": null, "overflow_x": null, "overflow_y": null, "padding": null, "right": null, "top": null, "visibility": null, "width": null } }, "a31d9848ee2a466fa3b8ee92ea83e9d2": { "model_module": "@jupyter-widgets/base", "model_module_version": "1.2.0", "model_name": "LayoutModel", "state": { "_model_module": "@jupyter-widgets/base", "_model_module_version": "1.2.0", "_model_name": "LayoutModel", "_view_count": null, "_view_module": "@jupyter-widgets/base", "_view_module_version": "1.2.0", "_view_name": "LayoutView", "align_content": null, "align_items": null, "align_self": null, "border": null, "bottom": null, "display": null, "flex": null, "flex_flow": null, "grid_area": null, "grid_auto_columns": null, "grid_auto_flow": null, "grid_auto_rows": null, "grid_column": null, "grid_gap": null, "grid_row": null, "grid_template_areas": null, "grid_template_columns": null, "grid_template_rows": null, "height": null, "justify_content": null, "justify_items": null, "left": null, "margin": null, "max_height": null, "max_width": null, "min_height": null, "min_width": null, "object_fit": null, "object_position": null, "order": null, "overflow": null, "overflow_x": null, "overflow_y": null, "padding": null, "right": null, "top": null, "visibility": null, "width": null } }, "d87aca4cf9444004b50a7d5176117880": { "model_module": "@jupyter-widgets/controls", "model_module_version": "1.5.0", "model_name": "SliderStyleModel", "state": { "_model_module": "@jupyter-widgets/controls", "_model_module_version": "1.5.0", "_model_name": "SliderStyleModel", "_view_count": null, "_view_module": "@jupyter-widgets/base", "_view_module_version": "1.2.0", "_view_name": "StyleView", "description_width": "", "handle_color": null } }, "e2d07d875ad34e97a02d11bcacd7c417": { "model_module": "@jupyter-widgets/base", "model_module_version": "1.2.0", "model_name": "LayoutModel", "state": { "_model_module": "@jupyter-widgets/base", "_model_module_version": "1.2.0", "_model_name": "LayoutModel", "_view_count": null, "_view_module": "@jupyter-widgets/base", "_view_module_version": "1.2.0", "_view_name": "LayoutView", "align_content": null, "align_items": null, "align_self": null, "border": null, "bottom": null, "display": null, "flex": null, "flex_flow": null, "grid_area": null, "grid_auto_columns": null, "grid_auto_flow": null, "grid_auto_rows": null, "grid_column": null, "grid_gap": null, "grid_row": null, "grid_template_areas": null, "grid_template_columns": null, "grid_template_rows": null, "height": null, "justify_content": null, "justify_items": null, "left": null, "margin": null, "max_height": null, "max_width": null, "min_height": null, "min_width": null, "object_fit": null, "object_position": null, "order": null, "overflow": null, "overflow_x": null, "overflow_y": null, "padding": null, "right": null, "top": null, "visibility": null, "width": null } }, "ed302fda8d554243a2ffa1c0decee4af": { "model_module": "@jupyter-widgets/controls", "model_module_version": "1.5.0", "model_name": "IntSliderModel", "state": { "_dom_classes": [], "_model_module": "@jupyter-widgets/controls", "_model_module_version": "1.5.0", "_model_name": "IntSliderModel", "_view_count": null, "_view_module": "@jupyter-widgets/controls", "_view_module_version": "1.5.0", "_view_name": "IntSliderView", "continuous_update": true, "description": "n", "description_tooltip": null, "disabled": false, "layout": "IPY_MODEL_3382016e18834c45b5b14ea3e427bdf1", "max": 10, "min": 0, "orientation": "horizontal", "readout": true, "readout_format": "d", "step": 1, "style": "IPY_MODEL_d87aca4cf9444004b50a7d5176117880", "value": 0 } } }, "version_major": 2, "version_minor": 0 } } }, "nbformat": 4, "nbformat_minor": 5 }